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ABSTRACT 

The   Harris    expansion    method   is    applied   to    the   elastic    s-wave 
scattering   of    low   energy    electrons    from   hydrogen    atoms    and   singly 
ionized  helium   atoms.      The    trial   wave    functions    are  Hylleraas 
functions    of   22,    34,    50   and   70  parameters.       It    is    found   that   rea- 
sonably   accurate   values    of   e-H  phase    shifts    can    be   calculated  but 
that    e-He      phase   shifts    are   substantially    less    reliable.       It    is 
shown    that    the    Harris    method   gives    an    accurate    depiction    of   the 
location,    but    not    the  width,    of    the    scattering    resonances.       Singlet 
and   triplet    s-wave  phase    shifts    for    e-H  and   e-He      scattering   are 
compared   with    the    results    of   other    calculations   and   H     and   He   S 
state    energy    levels    including    the   auto-ionizing    levels    are  pre- 
sented and   compared   with    other    calculations    and   with    experiment. 
It    is    tentatively    concluded    that    the    Harris    method   does    not    work 
on    systems    whose    long   range  potential    is    of    the    Coulomb   type. 
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INTRODUCTION 


The  quantum  theory  of  scattering  has  had  a  long  and  fruitful 
history,  beginning  with  the  very  origins  of  quantum  theory.   Its 
extensive  application  to  atomic  collisions  is  indicated  by  the 
large  number  of  texts,  monographs  and  reviews  of  that  subject  pre- 
sently in  print  [l].   More  recently  the  wide  availability  of  high 
speed  electronic  computers  has  given  additional  impetus  to  the 
subject.   New  methods  have  been  devised  and  older  ones  refurbished 
to  take  advantage  of  the  computational  power  now  at  the  disposal 
of  virtually  every  investigator  in  the  field.   Parallel  to  the 
development  of  computer  technology  has  been  the  application  of 
powerful  variational  techniques  to  the  scattering  problem. 

Variational  methods  have  held  a  distinguished  place  in  mathe- 
matical physics  since  the  days  of  the  Bernoullis.   Perhaps  no 
technique  has  had  wider  application  and  thus  it  is  no  surprise  that 
it  has  proven  a  powerful  tool  in  scattering  theory  as  well.   Its 
first  application  to  atomic  scattering  by  Hulthen  in  1944  [2J, 
followed  shortly  by  Kohn  [3]  and  Lippman  and  Schwinger  L4J  has  led 
to  extensive  applications  to  all  areas  of  scattering  theory  L  5 J • 
One  difficulty  that  arose  with  the  Kohn  and  related  methods  how- 
ever ,  was  the  appearance  of  spurious  singularities  in  the  cal- 
culated phase  shifts.   While  not  fatal  to  the  method  they  were 
annoying  and  some  considerable  effort  has  been  expended  in  attempts 
to  explain  these  singularities  by  Schwartz  [6]  and  others  L7,8,9J. 
By  working  to  exploit  rather  than  eliminate  the  singularities, 


Harris  [lOj  developed  an  expansion  technique  that  has  yielded  very 
reliable  results  to  date.   While  this  scheme  involves  some  compu- 
tational simplifications  over  the  Kohn  method  it  too  is  not  without 
disadvantages,  which  are  dealt  with  in  some  detail  below. 

Formal  solutions  abound  to  the  scattering  problem,  ranging 
from  the  Born  approximation  and  partial  wave  analysis  to  the 
Green's  function  methods  developed  by  Schwinger  LllJ  ,  but  accurate 
solutions  to  "real"  problems  beyond  the  realm  of  potential  scat- 
tering are  relatively  few,  owing  principally  to  the  complexity  of 
the  problem  when  one  permits  the  scattering  center  to  have 
structure.   Even  the  simplest  of  "non-trivial"  scattering  problems, 
that  of  electrons  or  positrons  incident   elastically  on  ground 
state  hydrogen  atoms  is  already  a  "many  body"  problem  and  thus  not 
subject  to  solution  in  closed  form. 

It  is  the  purpose  of  this  thesis  to  investigate  the  application 
of  the  method  developed  by  Harris  to  the  problem  of  the  elastic 
scattering  of  low  energy  electrons  from  hydrogen  atoms  and  from 
singly  ionized  helium.   Both  have  been  attempted  several  times 
previously  by  various  other  methods,  most  notably  in  the  case  of 
hydrogen  by  Schwartz,  whose  definitive  calculation  was  published 
in  196l  Tl2]  and  more  recently  by  Callaway  and  his  collaborators 
in  1969  and  1970  [l3,l4l,  and  by  Adelman  and  Reinhardt  [l5]  and 
Truhlar  and  Smith  Clo]  in  1972.   Because  of  the  more  complex  nature 
of  the  interaction  in  the  case  of  singly  ionized  helium  and  a 
paucity  of  experimental  data  for  comparison  these  calculations 
have  not  been  as  extensively  pursued.   The  earliest  attempt  was 
by  Bransden  and  Dalgarno  in  1953  Tl7]  and  the  most  fecent  by  Burke 
and  Taylor  in  1966  Cl8] . 
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With  the  exception  of  Schwartz'  work  roost  of  the  important 
calculations  of  these  problems  have  used  the  algebraic  close  cou- 
pling approximation  [l°]  wherein  the  trial  wave  function  is  made 
up  at  least  in  part  of  a  combination  of  the  first  few  bound  state 
orbitals  of  the  target  atom.   Such  a  scheme  has  earned  a  deserved 
reputation  for  accuracy  and  rapid  convergence,  however  from  a 
philosophical  point  of  view  it  suffers  because  its  application 
requires  some  a  pr  ior  i  knowledge  of  the  detailed  structure  of  the 
target  particle.   On  the  other  hand  Schwartz,  in  his  calculation, 
made  use  of  an  ever  increasing  set  of  those  functions  first  used 
by  Hylleraas  in  his  famous  calculation  of  the  ground  state  energy 
of  helium  L20j  .   This  method  makes  no  approximations  other  than 
those  required  by  the  use  of  a  finite  rather  than  an  infinite  basis 
space  in  which  to  work  . 

It  was  in  the  spirit  of  Schwartz'  work  that  the  present  project 
was  undertaken.   It  was  hoped  that  by  utilizing  the  Hylleraas  type 
functions  in  a  calculation  of  electr oh -helium  ion  scattering  that 
a  more  accurate  and  complete  compilation  of  s-wave  phase  shifts 
would  be  obtained.   As  will  be  seen  below  such  hopes  were  not 
realized,  however  it  is  hoped  that  the  reporting  of  such  results 


This  is  not  strictly  true.   Electron  capture  by  the  scattering 
center  is  usually  neglected  in  calculations  of  this  type,  and  is 
neglected  also  in  the  calculation  reported  in  this  thesis.   To 
incorporate  capture  requires  including  terms  in  the  potential  in- 
volving the  interaction  of  the  electrons  with  the  electromagnetic 
field  in  order  to  account  for  the  creation  of  the  necessary  photons 
and  is  properly  a  problem  in  quantum  electrodynamics.   The  reported 
cross  sections  for  capture  in  the  system  under  consideration  here 
are  small  compared  to  the  elastic  scattering  cross  sections  and 
usually  may  be  safely  neglected. 


as  were  obtained  will  provide  incentive  for  others  to  seek  an 
explanation  of  the  apparently  anomalous  results  for  helium  ion 
scattering  herein  reported. 

In  Section  II  a  brief  development  of  the  partial  wave  method 
of  scattering  theory  will  be  presented  including  the  modifications 
required  in  the  presence  of  the  long  range  Coulomb  force  followed 
by  a  description  of  the  Kohn  method.   Last  will  be  a  more  detailed 
description  of  the  Harris  method  and  its  application  to  elastic 
electron  scattering  by  single-electron  atoms  and  ions.   In  Section 
III  the  results  of  the  numerical  calculations  are  presented  for 
both  atomic  hydrogen  and  singly  ionized  helium  along  with  a  dis- 
cussion of  these  results  and  comparison  with  those  of  other  workers 
An  important  by-product  of  the  Harris  method  is  the  calculation  of 
the  bound  state  energy  levels  of  the  corresponding  two-electron 
system.   Results  of  these  calculations  for  H   and  He  are  also 
presented  in  Section  III.   Most  of  the  detailed  calculations  as 
well  as  discussion  of  the  numerical  methods  employed  are  deferred 
to  the  appendices. 
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II.   EXPANSION  TECHNIQUES  IN  SPATTERING  THEORY 

In  this  section  will  be  presented  a  brief  outline  of  the  method 
of  partial  wave  analysis  [21]  for  the  case  of  elastic  atomic  scat- 
tering.  Part  B  will  briefly  present  the  Kohn  method  and  in  part  C 

the  Harris  method  will  be  developed  in  some  detail  along  with  its 

+ 
application  to  elastic  He   scattering  of  electrons.   The  system  of 

units  used  throughout  will  be  those  in  which  m  (electron  mass), 

ti     (Planck's  constant)  and  e  (electronic  charge)  all  equal  unity. 

In  such  a  system  the  unit  of  energy  is  equal  to  twice  the  ground 

state  energy  of  hydrogen  (i.e.,  27.2  eV),  and  momenta  are  written 

2 

in  terms  of  the  wave  number,  k,  such  that  E  =  k  /2 .   It  will  be 

assumed  that  the  nuclear  mass,  be  it  hydrogen  or  helium,  is  infi- 
nitely large  compared  to  that  of  the  electrons. 

A.   PARTIAL  WAVES  AND  PHASE  SHIFTS 

Let  a  beam  of  particles  of  momentum  k  fall  upon  a  spherically 
symmetric  potential  V(r).   Taking  the  scattering  center  as  the 
coordinate  origin  and  the  direction  of  k  as  the  angular  axis  the 
wave  function  of  the  outgoing  beam  far  from  the  coordinate  center 
will  be  proportional  to  the  sum  of  the  wave  function  of  the  inci- 
dent  beam,  represented  by  the  plane  wave   exp(ik-r)  and  a  modified 
spherical  wave,  f (9) exp ( ikr )/r ,  representing  the  scattered  portion 
of  the  incident  beam.   For  proof  that  solutions  to  the  Schroedinger 

equation  of  the  form 

— »  — * 
Y(*)^eik'r  ♦  f(0)eikVr  (1) 
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exist,  the  reader  is  directed  to  reference  21  (p.  220  ff).   The 
differential  scattering  cross -sect ion  is  defined  to  be 

Number  of  particles  per  unit  time  scattered  into 

.  .    -    the  solid  angle  element  d-a  at  .a 

*  '     A  ~  Total  Incident  Flux  '  ' 

where  the  particle  flux  is  defined  to  be 

q^F"  (Y*VY  -  YVY*)  .  (3) 

Putting  (3)  into  (2)  gives 

-*         -.        -»      A       2 
cp    -ds        cp    *r    r 

CT(9)da   =  — =  ~ da  (4) 

v      '  CP  CD  v      ' 

o  o 

-  A 

where   ds    is    the    element    of   surface   area   at    the   detector    and   r    is    a 

unit    vector    in    the   direction    of   the   scattered   flux.      cp      and  cp 

o  s 

represent  respectively  the  flux  arising  from  the  first  and  second 
terms  of  (1),  respectively.   Because  of  the  scalar  product  in  (4) 
only  the  radial  part  of  the  scattered  wave  function  contributes 
to  the  scattered  flux.   Evaluating  cp   and  cp   from  (3)  and  applying 
them  to  (4)  leads  immediately  to 

a(9)drL  =  |f (9)l2  da  .  (5) 

It  remains  then  to  exhibit  a  means  of  calculating  f(9)  to 
complete  the  connection  between  theory  and  experiment.   Such  a 
means  is  the  method  of  partial  wave  analysis.   The  most  general 
solution  to  the  Schroedinger  equation  in  the  case  of  axial  sym- 
metry is 

Vr)   o 
1=0      *         r     Z 
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where   Y*(Q)    is    the    spherical    harmonic    of    order    zero.      The    index  i 
represents    the  angular    momentum    quantum   number.      Each    radial    function 
satisfies    an    equation    of   the   form 

If  V(r)    goes    to   zero   fast    enough   at    large   r    there   will  be  a    region 
outside    of   which    its    effects    can   be   neglected.      Since   u.(r)    is    a 
function    of  r    only,    equation    (7)    can   be   rewritten 

u"    +    [k2-2V(r)-   £(£+l)/r2]    u    =   0 

and    since    the   last    two   terms    tend   to  zero   for    large    r    it    is    rea- 
sonable   to  assume    that   any    solution    u   tends    to   the   form  A*sin(kr+e) 
for    large   r.       If    this    be   so,    u   must    have    the   form   f(r)e  where 

f(r)    is    approximately    constant    for    large   r,    and    (7)    reduces    to 

f"    +    2ikf    -    [2V(r)+   £(£+l)/r2]    f   =   0 

and    for    large    r,     if    f    is    nearly    constant    f  «    kf      and   this    equation 
can    be  approximately    integrated,    giving 

2ik&i  f    =    U2V(r)+   £(£+l)/r2]dr 

and    the   right    hand   side    tends    to   a    constant    for    large   r    only    if 
V(r)    goes    to    zero   faster    than    1/r     (ref.    22,    p.    23). 

In    the    region    where  V(r)    can    be   neglected   equation    (7)    reduces 
to   the   spherical    Bessel    equation   whose   solutions    are 

— =   A^kr)*    B^n^kr)  (8) 
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where   j„(kr)    and   n . ( kr ')    are    the   spherical    Bessel   and  Neumann 
functions,    respectively    L23J    whose   asymptotic    forms    are 

J£(k*)^rZ>  sin(kr-£TT/2)/kr 

nje(kr)^zr5  cos(kr-je  n/2)/kr   . 

2 
Since  n. (kr)  is  not  well  behaved  near  the  origin  ,  if  V(r)  =  0 

everywhere  then  B.  must  be  chosen  to  be  zero  in  order  for  the 

solution  u,  to  be  valid.   Thus  it  can  be  inferred  that  the  value 

of  B.  compared  with  A-  will  be  a  measure  of  the  strength  of  the 

potential  and  hence  of  the  scattering.   To  make  this  explicit  (8) 

can  be  rewritten,  using  the  asymptotic  forms  above 

A„  r         *_    B, 

k 

h 

and  letting  B-/A.=  -  tan6. 


u^r)  —  T   [sxn(kt-  -)-   T    cos(kr-  -)  j 


Vr>F^CXsin(kr'  -j-  +  6^)  (9) 

so  that  equation  (6)  becomes  for  large  r,  a  plane  wave  shifted  in 
phase  from  that  of  an  unperturbed  (V=0)  plane  wave.  C.    can  be 
taken  equal  to  one  and  the  normalization  absorbed  by  the  expansion 
coefficients  in  equation  (6). 

The  first  term  on  the  right  hand  side  of  (1)  can  be  expanded 
as  follows: 

00 

eikrcos9  =  s=o(2i+l)i£jjt(kr)Pi(cos©)  (10) 


2 

The  Neumann  function  goes  as 


r      r 
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and  letting  j.    assume  its  asymptotic  form,  (6)  and  (1)  can  be 
equated  and  the  result  solved  for  f(0),  giving,  on  choosing  the  a, 
so  as  to  assure  an  outgoing  wave, 

i6 


f(9)  :  J7  £  (2l+l)e      ^sin6  P.(cos9) 
K  1=0  l   l 


(11) 


which  is  the  usual  form  given  for  the  scattering  amplitude  in  terms 
of  the  phase  shift. 

Deferring  the  calculation  of  6.  to  the  next  sections,  the  sub- 
ject of  how  the  above  formalism  must  be  modified  in  the  presence 
of  a  long  range  coulomb  potential  will  now  be  considered.   Recall 
that  the  classical  solution  for  the  motion  of  a  particle  in  an 
inverse  square  central  force  field  is  a  hyperbola.   According  to 
an  argument  due  to  Gordon  (ref.  22,  p.  54),  if  one  considers  the 
family  of  classical  hyperbolic  trajectories  with  one  asymptote 
originating  at  the  left  limit  of  the  z    axis,  the  surface  perpen- 
dicular to  these  hyperbolae  goes  as 

i 
Z+  a  \  k(r-z)  =  const 

for  large  r.   Thus  it  is  as  if  the  incident  wave  were  initially 
distorted  by  the  presence  of  the  scattering  center  at  infinity  and 
its  expected  form  would  then  be 
ik{z+  a  07!  k(r-z) } 

A  similar  argument  with  regard  to  the  scattered  wave  leads  to  the 
form  for  the  asymptotic  wave  function  in  the  presence  of  a  coulomb 

field 

(12) 


...  .  .      ikz  +  iafo  k(r-z)  -,n\    ikr 
Y(r)— -J  e  v    '+   f(9)e  . 


-  ia.1-.  2kr 
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The  coefficient  a  will  be  defined  below.   In  solving  equation  (7) 
in  the  presence  of  the  Coulomb  potential  one  typically  converts  to 
parabolic  cylindrical  coordinates  whereupon  the  Schroedinger 
equation  becomes  a  version  of  the  hyper geometr ic  equation  whose 
solutions  are  the  hyper geometr ic  functions  (ref.  22,  p.  57  f f ) . 
Using  the  asyrapototic  forms  of  the  solution  functions  as  before, 
equation  (9)  becomes,  in  the  case  of  the  Coulomb  potential 

in 
u£(x)Tm  C£sin(kr-  -j   +  a  Bn   2kr+  r\^   6^)  (13) 

where    6.    is    as    defined  above   and 

77 £    =    arg  [T(i+1    -    ia)}  (14) 

is  the  phase  shift  at  infinity  due  to  the  Coulomb  potential. 
Equation  (11)  must  now  be  modified  to  account  for  this  additional 
phase  shift  (ref.  22,  p.  65  ff). 

The  coefficient  a  in  (12),  (13)  and  (14)  is  proportional  to 
the  charge  of  the  scattering  center.   However,  since  the  problem 
at  hand  is  one  in  which  the  pure  Coulomb  field  exists  only  at  long 
range,  this  charge  is  the  net  charge  of  the  ion,  in  this  case  (Z-l) 
where  Z  is  the  atomic  number  of  the  scattering  center,  and  in  the 
system  of  units  in  use  a  can  be  written 

a  =  (Z-l)/k  (15) 

and  now  5   can  be  interpreted  as  the  phase  shift  due  to  the  de- 
parture from  pure  coulomb  scattering  at  close  range.   It  is  this 
quantity  which  is  of  physical  interest. 

To  summarize,  it  has  been  shown  how  the  asymptotic  form  of  the 
wave  function  (equations  1  or  12)  can  be  calculated  by  solving  the 
Schroedinger  equation  for  a  series  of  linearly  independent  functions 
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(equations  9  or  13),  in  other  words,  that  the  problem  can  be  reduced 
to  solving  a  series  of  "partial  wave"  equations  for  the  various 
values  of  the  angular  momentum  quantum  number  i .      A  simple  kine- 
matic argument  indicates  that  as  the  energy  of  the  incident  beam 
is  reduced  the  higher  order  terms  in  the  expansion  (6)  become  less 
important,  enabling  a  fairly  good  representation  of  f(9)  to  be 
constructed  using  only  the  first  few  terms. 

B.   THE  KOHN  VARIATIONAL  METHOD 

Before  taking  up  the  Harris  method  it  will  be  worth  while  to 
briefly  examine  the  Kohn  method  since  it  is  the  most  widely  used 
variational  method  in  low  energy  atomic  scattering  and  its  dif- 
ficulties were  the  motivation  for  Harris'  development.   These  dif- 
ficulties have  been  examined  elsewhere  [6-9J  and  will  only  be 
mentioned  here. 

The  Kohn  method  provides  an  estimate  of  the  phase  shift  which 
is  "second  order  accurate"  in  the  error,  i.e.,  the  error  term  is  at 
least  second  order  small.   In  developing  the  Kohn  method  one  assumes 
a  trial  function  of  the  form 

(V    =   sin  kr  +  t  cos  kr  +  Y  (16) 

where    t    is    an    estimate    of    the    tangent    of   the  phase   shift    and  Y    is 

normally   a    linear    combination    of    some   basis    set   X-    such    that 

N 
Y    =  S        a.\-    -ZT2>°  (17) 

so   that  cp    is    now   a    function   of   N+l   parameters—the   N   values    of  a., 

and   t.      Let   ^  be    the   true    wave   function    whose   asymptotic    behavior    is 

^    — -J  sin    kr+    tan5    cos    kr 
r-  ro 

¥(0)=  cp(0)    =   0 
17 


and  define  the  error  wave  function 

e  =  cp  -  ^.  (18) 

Now  form  the  functional 

F  =  \  cp(H-E)cp  dr  (19) 

which,  using  the  fact  that  (H-E)^©,  can  be  transformed  into 

F  =  V  cp(H-E)e  dr  . 

Integrating    twice   by   parts    yields 

p  =    _  -(cpve    -   eVcp)  +  V    e(H-E)cp    dr 

and  substituting  (l8)  for  £  in  the  first  two  terras  and  for  cp  under 
the  integral  gives 

F  =  \  Cp(H-E)<p  dr  =  -(t-tan5)+  \  G  (ll-E)G  dr 


or 


|  tan5  =  |  t-  F  +  \  6(H-E)e  dr  (20) 

which    identity    is    due    to  Kato    V2l\~\    and   exhibits    the    nature    of    the 
error    terra.      Now    the   Kohn   prescription    takes    as    its    functional 

I=|  tanS    =   |  t    -    \    Cp(H-E)cp    dr*  (21) 

and  seeks  its  stationary  value  subject  to  the  conditions 

bl    &i 


bt    ba. 

x 


=  0,      1=1, . . . . ,N. 
The  value  of  I  thus  found  is  the  estimate  of  tan6  to  second  order 


C.   THE  HARRIS  EXPANSION  METHOD 
1 .   Formal  Solution 

It  can  be  shown  that  the  development  of  expressions 
for  the  phase  shift  lose  no  generality  when  the  system,  is  restricted 
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to    one    of    s-waves    (-£=0)    scattered   from   a    central   potential.      There- 
fore   this    restriction   will    be    imposed   on    the  present    derivation   of 
Harris'    method   Lio]    and   the    validity    of    its    extension    to  atomic 
systems    will    be   assumed. 

With    the   above    considerations    in   mind,    assume   a    system 
describable   by   a  potential   V(r)    and  a   Harailtonian   H.       Construct   a 
trial    basis    space   using    a    set    of   N    linearly    independent    functions 
iX-^'       I"    this    restricted   space    the   Harailtonian    can    be   represented 
by   a   matrix   whose   elements    are 

Hij    =  <XjlIhIx  >  (22) 

while   the    inner    products    of   the    basis    vectors    are 

Lij    =  <Xi'XJ>       '  (23) 

Then    the    eigenvalue   equation 


(H    -  \L)C    =    0 

where  the  elements  of  H  and  L  are  as  given  in  (22)  and  (23),  can 
be  solved  for  the  N  eigenvalues  A.   and  corresponding  eigenvectors 
C  ,  (p  =  l  ,  .  .  ,N)  .   From  the  eigenvectors  a  new  basis  set  is  con- 
structed (which  spans  the  same  space  as  the  IX^) 


] 


N 


V  =  \,1ciP]xi>  ■  (25) 


Now  if  the  true  wave  function  is  represented  by 

^  =  S+tC  +  $ 

where  ^(0)=  0 

and  T 

^(r  )  -n'sin  kr  +  tan6  cos  kr  =  S  +  tC 


19 


(where   now    the    identification    of    the    terras    S,    t    and  C    is    clear) 
and    if   in   addition    \§>   may    be   reasonably    well    represented  at    some 
energy   E  by    \^>> ,    a    linear    combination   of    the    lcPn>> 

N 

k>  =  r     b  |  co  > 

then  the  approximate  wave  function  may  nearly  satisfy  the 
Schroedinger  equation.   If  this  is  to  be  the  case,  then  in  the 
space  spanned  by  the  l9n>  the  vector 

(H-E) |s+  tC+  5>  (26) 

must  be  the  null  vector,  which  is,  of  course,  merely  the  statement 
of  Schr oedinger ' s  equation  in  the  restricted  space  of  the  trial 
functions.   This  condition  may  be  satisfied  if  it  is  required  that 
(26)  have  no  component  in  the  space  spanned  by  the  | tp-5* >  that  is, 
that 

<cp  |(H-E)|S+  tC+  l>    =  0,      p=l,...,N  .  (27) 

Since  the  cp  >    are  themselves  linear  combinations  of  the  original 
0 

basis  set  IX-^j  "the  condition  (27)  is  equivalent  to  requiring  that 

(26)  have  no  component  in  the  space  spanned  by  the  |X.>«   Bearing 

in  mind  that  the  solution  sought  involves  finding  t,  an  estimate  of 

tan6  ,    note    that    if  <(D    |h-e|§>   and  <cp    |h-e|s+    tC>   are   separately 

p  p  ^ 

equal  to  zero  then  (27)  is  satisfied  and  t  can  be  immediately  found, 
Note  further  that  if  E  =  X     ,  the  eigenvalue  found  in  (24)  corre- 
sponding to  the  vector  |cPn>,  that  then 

«Ppl(H-E)|§>  =  0  (28) 
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for 

N  N  N 


«pJ(H-E)|§>    =   E         b      S       C        S      <X     |(H4  D)|X,>C 
P  w=l    w  j=l   jCc  i=i      x  p        J 


ip 


and   the    third   sum   in   this    expression    is    merely    (24)    in    component 

notation,    whence 

N 

£       <X-I  (H-X     )|X  •>    c         =    0,  j=l,...,N 

i=l  ^  J  H 

irrespective    of    the   value    of   CO.       Hence    (28)    is    satisfied  and    it 
follows    immediately    from    (27)    that 

<co    I  (H-E    )|s> 
ta"6    =    t    =    "    <V(H^)|C>  <29' 

and  the  problem  is  formally  solved.   Finally  as  the  size  of  the 
basis  set  |X->  is  increased  (i.e.,  approaches  a  complete  set  of 
quadrat ically  integrable  functions)  |  §>  should  more  closely 
approach  | I  > and  t  should  converge  to  the  correct  result. 

Kolker  [25]  has  shown  that  the  Harris  method  can  be  for- 
mulated as  a  variational  principle  and  that  under  certain  con- 
ditions it  can  be  combined  with  the  Kohn  method  to  give  a  minimum 
principle.   The  advantage  of  the  Harris  method  is  principally  com- 
putational, in  that  no  integrals  of  the  form  <S|  (H-E)|C>  need  be 
evaluated.   However  the  contribution  of  these  integrals  is  needed 
to  provide  the  second  order  correction  to  the  estimate  of  tan6  in 
variational  methods  dependent  upon  the  Kato  identity  (20),  so  it 
can  be  seen  that  the  Harris  method  is  necessarily  less  accurate 
than  the  Kohn  method,  as  was  found  by  Nesbet  [8].   It  remains  true, 
however,  that  the  correction  term  can  be  made  arbitrarily  small  by 
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increasing   the   size   of    the    trial    basis    set,    for    as    this    set 
approaches    completeness,    the    error    roust   approach    zero. 
2 .       Application    to   Atomic    Scattering 

In   applying   the   Harris    formalism    to   the  problem    of   atomic 
scattering    the   Hamiltonian    is    written   as    follows 

12        12       7        7  1 

2122rrr  KO    ' 

1         2         12 

where   the    subscripts    refer    to   the    individual    electrons,    Z    is    the 
atomic   number    of   the   scattering   atom  and   r         =    |r       -    i     |     is    the 
inter -electr on    distance.      The   zero   of   total    energy    in    this    system 
occurs    with    all    three   particles    at    rest    and    separated    from   one 
another    by    an    infinite    distance.      The   trial    basis    set    is    taken   as 
the   Hylleraas-type    functions    [20] 

(    n    i  I    n\    m  -  — (  r  ,+r „ )  ,„,. 

Xi    =    Vir2    -  tir2>'ri2    e      2\    !       2y  (31) 

where   n,    t    and   ra  are    integers,    chosen    subject    to    the    constraint 

n  +  I   +  m^M,  (M=l  ,2,  .  .  .  )  (32) 

and  a  is  a  variable  parameter. 

This  form  of  the  wave  function  allows  for  electron  exchange, 
i.e.,  the  interchange  of  roles  between  bound  and  free  during  the 
scattering  process,  and  enables  the  wave  function  to  be  properly 
symmetrized  when  the  electron  spins  are  anti-parallel  (spin  =  0, 
or  "singlet"  state)  or  anti-symmetr ized  when  the  spins  are  parallel 
(spin  =  1 ,  or  "triplet"  state)  thus  satisfying  the  requirements  of 
the  Pauli  principle  [26,27].   This  is  the  only  way  the  spin  enters 
the  problem  and  so  once  the  symmetry  of  the  wave  function  is  prop- 
erly chosen  no  further  concern  need  be  taken  with  the  spin 
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coordinates  of  the  electrons.   Inclusion  of  the  r    term  in  the 
potential  and  in  the  wave  function  takes  account  of  the  induced 
polarization  of  the  atom  arising  from  the  repulsion  felt  by  the 
bound  electron  as  the  incident  electron  approaches.   Finally,  the 
exponential  term  forces  the  functions  to  zero  as  r   or  r   increases 
so  that  all  boundary  conditions  required  for  these  functions  in  the 
last  section  are  met.   As  will  be  seen  below  the  parameter  d.    is 
varied  by  inspection  so  as  to  optimize  the  results  obtained. 

The  number  of  functions  in  the  basis  set  is  found  by  taking 
all  possible  combinations  of  n ,  t,    m  subject  to  the  constraint  (32) 
for  a  given  M,  and  eliminating  all  terms  which  duplicate  the 
function  or  require  it  to  be  identically  zero  irrespective  of  the 
choice  of  sign.   The  results  for  1^  m  <.   8  are  shown  in  Table  I. 


M  No.  terms  in  basis  set 

Singlet  Triplet 

1  3  1 

2  7  3 

3  13  7 

4  22  13 

5  34  22 

6  50  34 

7  70  50 

8  95  70 


Table  I.   Basis  set  size  for  various  M. 

The  choice  of  coordinate  system  is  of  importance  in  sim- 
plifying the  computational  details.   For  the  present  computation 
the  optimum  choice  seemed  to  be  the  one  in  which  the  principle 
coordinates  are  the  radius  vectors  of  the  two  electrons  and  the 
angle  formed  by  the  radius  vectors.   The  six  coordinates  governing 
motion  of  the  center  of  mass  and  orientation  of  the  system  in  space 
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are    cyclic,    and   the   Harailtonian   becomes    L28J 


».  2      6  b  2  __ft_     /  J_        _1 


sinO 


12    bO 


12 


2  Z     )       1 

r  r  r 

1  2  12 


(33) 


where   r^    =   (r*    +    r2    .   2x^ 

the   volume   element    for    integrals    in    this   system    is 


r       cos9      J      by    the   law    of   cosines,    and 


d7ldr2   =   8n2«Jdr1r|dr2sin012dB12    . 


(34) 


The  advantage  of  this  system  is  that  all  of  the  integrals  arising 

from  forming  the  matrices  H  and  L  (equations  22-24)  can  be  done 

exactly,  while  those  arising  from  <cp  |  ( H-E)  |  S  +  tC  >  can  be  done 

exactly  in  the  case  of  hydrogen  and  reduce  to  single  integrals 

which  then  can  be  done  numerically  in  the  case  of  helium  ions. 

Because  the  electrons  are  identical  and  thus  formally 

indistinguishable  it  is  not  necessary  to  integrate  over  the  entire 

first  quadrant  of  the  r  ,r   plane.   It  is  sufficient  to  cover  only 

the  region  r  ^r   [29].   To  see  why  this  is  true  one  must  examine 

the  integrals  generated  in  formino  H.  .  and  L.  .  usinq  (31)  and  (33) 

xj  xj  v   '      v 

Carrying  out  the  indicated  operations  gives  rise  to  52  terms  for 

Each  term  consists  of  an 
integral  of  the  form 


each  H.  .  and  four  terms  for  each  L.  . 


A(V,X  ,u)     =   \  dr    dr 


"   *(ri+r2>       V-2   X-2   H-2 
2    e  ri       r2      ri2 


(35) 


or 


B(v,X  ,H)     =   W^ 


-a(ri+r2)       V-2X-2U-2 
r2    e  ri       r2      ri2      C°S912    ' 


(36) 
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The  complete  expressions  for  H.  .  and  L.  .  are  given  in  Appendix  A. 

Examination  of  the  terms  comprising  H.  .  and  L.  .  indicates 

that  every  term  involving  A(V,X,0)  and  B(v,A>0)  always  has  a 

coefficient  of  zero.   Since  the  angular  integration  of  B(V,A>2) 

vanishes,  B(v,A,2)  =  0,  so  using  the  recursion  relations  proved 

in  Appendix  B  it  is  sufficient  to  explicity  evaluate  only  the 

A(v,X,2)  in  order  to  find  all  the  terms  required  to  generate  the 

H. ,  and  L  .  .. 
ll       ij 

Setting  a.    =    2    in    (35)    and   performing    the  angular    integration 
over    the    range   0    to  tt  ,    gives,    in    terms    of   the    integration    region 
outlined    in    Figure   1    below   and   the    volume    element    (34) 


P   °°  -   ar  ^1 

A(V,X  ,2)    =    16  tt2  \         r^    e  dr      \         r 

J0  J0 


2    G  *2    '  (37) 


Inspection  of  H.  .  shows  that  for  each  term  of  the  form  (37)  there 
is  another  of  the  form 


A(X  ,V,2)     =    16  rr 


.  °°  %       -  ar,r.    1  -  Q!r 

driri    e  \         dl2r2    6 

0  J0 


(38) 


which  is  the  mathematical  equivalent  of  integrating  (37)  in  the 
region  r   ^  r   as  can  be  seen  immediately  by  making  the  variable 
change  r  ~   r   in  (3S).   Thus  as  a  result  of  the  exchange  of 
particles  possible  because  the  electrons  are  identical,  the  in- 
tegrals (37)  are  sufficient  to  cover  the  entire  quadrant.   This 
argument  is  extended  below  to  the  integrals  arising  in 
<cp  I  (H-E)  |s+tC>  . 

Equation  (37)  may  be  integrated  directly,  giving 


,   ,  „x    16  tA 

A(V,X,2)  = — -— 

V+X+2 


v.   lj   (v+X-iH, 


2  ^=0  2V+X"1(X-i)I  . 


(39) 
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Figure  1.   Integration  Region  for  A(v,A>2)« 


2 
The  factor  l6  n  appears  in  every  integral;  it  can  be  dropped  from 

the  calculation  since  it  will  merely  give  rise  to  the  same  constant 
factor  in  both  H  and  L  and  will  not  affect  the  solution  of  the 
problem.   Because  of  the  way  the  recursion  relations  are  structured 
the  factor  \/0L   always  occurs  in  the  form  1/0!       and  thus  can  be 
factored  out  of  the  integration  and  incorporated  into  the  coeffi- 
cient of  each  term  in  the  expression  for  H.  .  and  L.  ..   Since  this 
is  the  only  way  the  parameter  CL   enters  these  integrations,  evalu- 
ating all  the  integrals  need  only  be  done  once  at  the  beginning  of 
the  calculation  for  a  given  basis  and  the  values  thus  found  can  be 
used  over  again  as  OL    is  varied. 

Unfortunately  a  similar  simplification  is  not  available  with 
the  integrals  from  <  cp  |  (H-E)  |  S+tC  > .  However  as  will  be  seen  only 
a  relatively  few  integrals  must  be  done  numerically  for  each  value 
of  0L   and  k  and  so  the  time  required  to  complete  the  numerical 
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integration  (which  dominated  the  time  required  to  evaluate  the 
helium  ion  phase  shifts  in  all  but  the  largest  basis  sets)  was  kept 
within  reasonable  bounds.   Since  the  present  calculation  is  con- 
cerned only  with  the  s  partial  wave  (-£=0)  further  discussion  will 
assume  i    has  been  set  equal  to  zero  in  all  equations  of  the  form 
(9),  (13)  or  (14).   Furthermore,  since  the  proper  form  for  hydrogen 
results  if  (13)  is  taken  as  the  asymptotic  form  of  the  wave  function 
and  Z  is  set  equal  to  one  in  (15),  the  remainder  of  the  development 
will  be  concerned  exclusively  with  establishing  the  Harris  method 
in  the  presence  of  a  long  range  Coulomb  field. 

Bearing  these  comments  in  mind,  note  that  as  r  goes  to 
zero  (13)  does  not  behave  as  well  as  one  would  like.   Because  of 
the  presence  of  the  logarithmic  term  in  the  argument  of  the  sine, 
as  r  approaches  zero  the  interval  between  zeroes  of  the  function 
approaches  zero  and  u(r)  does  not  approach  a  definite  limit.   To 
force  u(r)  to  zero  a  factor  is  normally  chosen  which  goes  to  zero 
fast  enough  as  r  goes  to  zero  to  override  the  effects  of  any 
singularity  which  may  occur,  and  goes  to  unity  as  r  gets  large  in 

order  not  to  affect  the  asymptotic  behavior  of  the  function.   For 

3 
the  present  calculation  the  factor  was  chosen  to  be  (1 -exp( -Wr /2) )  . 

This  choice  was  made  to  insure  that  u(r)  would  go  to  zero  fast 

enough  that  its  oscillatory  behavior  near  the  origin  would  not  affect 

the  accuracy  of  the  numerical  integrations  involving  the  asymptotic 

part  of  the  wave  functions. 

With  these  preliminaries  aside,  the  form  chosen  for  the 

asymptotic  part  of  the  wave  function  may  now  be  presented: 
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a 
s+tc  =  G 


(l-e      2      2)3[sin(kl'2+   ^JT^   2kr2+  ^ 


tan6    cosfkr    +  -^—£7:  2kr    +  77^ J/r 


-Zr,         --, 

+   e 


(l-e      2      1)3[sin(krl+  -~  Urn.  2kr1+  T|) 


where 


+   tan6    cos(kr    +  ~i  2m  2kr    +  r?)j/r  (40) 

7J    =    arg  {r(l-i   ^-  )}     .  (4l) 


Evaluation    of   the    complex   Gamma    functions    is    discussed    in   Appendix 

D.       Normalisation    factors    in    (40)    have   been   suppressed   because   they 

do   not    affect    the    calculation,    cancelling   when    t    is    evaluated. 

Following   the  prescription    outlined    in    section    1,    the 

estimate    of    tan5    is    evaluated   from   theratio   of    the    integrals 

<  ®D  l(H-E)  I  S>  and  <  qp      |  (  H-E)  |  C>     at    the    energy    E      found   by    solving 
r  r  r 

(24).       Since    the   functions    (H-E)|X->    have   been   previously    evalu- 
ated   in    solving    (22)    and    (23),    considerable    simplification    is 
achieved   by    taking   advantage    of   the    fact    that    H    is    a    hermitian 
operator    and    solving    instead   the    equivalent    integrals 

N 


U  cip<s,<H-y|xi: 


(42) 


and 

N 

Z        c       <c|(h-e   )|x.>  (43) 

i=l      XP  P         x 

whence  t  is  found  from  the  negative  ratio  of  (42)  to  (43).   Note 
again  that  6  measures  only  the  departure  from  pure  Coulomb 
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scattering   arising    in    the  region    close    to   the   target    where   the 
three  particles    in    the   system   must   be    treated  explicity. 

Carrying    out    the    indicated  algebra    one  finds    that 
<S|(H-E)!\>    has    52    terms    constructed    of    integrals    of    the   form 

,      ,        .       C   V-2   X-2   u-2       "(Z+2,rL       "2I2(1       ~2r2) 
A(v,X,Li)=y         r         r  e  e  M-e  /         x 


Z-] 


sin(kr    +  ~—  m  2kr    +Tj)  dr    dr 


(44) 


A2(v,X,u)=   ^ 


-2  X-2    u-2 
r2       ri2       e 


a 
-(Z+o)r. 


"  2  ri 
e 


1-e 


2 


^) 


Z-l 


sin( kr    +   — —  #rc   2kr    +  77)  dr  ,  dr  ^ 
v       1         k  1      ' '        1      2 


(45) 


^.j,     in_fv-2X-2U-2       -(Z+^ri"2r2       (  -2r2) 

B   ( v,X  ,[i)=   \r         r  „      r ,  „      e  e  \l-e  /         x 


r2      ri2       e 


Z-] 


sin(kr2+  —  0n  2kr2+  T))cos9      dr    dr  (46) 


B2(V 


'x^>=Sri 


"(Z+    2^2    "2ri 
r2       ri2       e 


-2    X-2    H-2 


1-c 


a 

2 


-   \3 
/      x 


Z-l 


Sln(kri+   ~~k~  ^  2kr;L+  r7)cos912dr1dr2       (47) 


Four    similar    integrals    with    cos(kr    +    .  ...)    replacing   the   sine    terms 
in    (44)    -    (47)    arise    in    evaluating  <    c|(H-E)|x>     and  are    designated 
A^,    A^,    B    ,    and    B, ,    respectively.      Note    that   all    these    integrals 
are    of    the   same    form  as    those    in   Appendix    B  Section    1   and    hence   the 
same    recursion    relations    may   be   applied. 
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As  with  the  integrals  of  equation  (34)  the  evaluation  of 
(44)  -  (47)  starts  with  the  case  U-  =  2  and  as  before,  the  angular 
integration  yields  B  (v,\,2)  =  B  (v,X,2)  =  0  and  introduces  a 
factor  2  in  the  coefficient  of  A,(v,X,2)  and  A  (v,X,2).   Also  as 
with  the  A(v,X,2),  terras  of  the  form  A.(V,\,2)  and  A„(v,X,2) 
always  exist  together  and  thus  satisfy  the  requirements  that  allow 
the  integration  to  proceed  over  r   ^  r   only,  as  illustrated  in 
Figure  1 . 

Starting  first  with  A  (v,X,2),  after  angular  integration 

Al(v,X,2)=  16^2  I   "dV^  ^'^"A''   ^l-"^2)3  x 

sin(kr2+  ~^   2kr2+  ^  ' 

Now   referring    to   Figure   2,    note    that    the    order    of    integration    can 

be   reversed  by    choosing    the   elements    of   area    dr    dr      as    indicated 

while    the   same    region    is    covered.      With    this    change   A,(V,X,2) 

becomes 

a        ,  a 


A;L(v,X,2)    -    l6^r2   Jj    dr2r^e      2      2Vl-e      2    ' 2)         x 

sin(kr    +  ~-  0n  2kr2+  r})\    dr    r^e-*Z+2  )ri 


r2 


and  now  the  integral  over  r   can  be  done  giving 


«   \3 


A  (v,X,2)  =  ^r-r—  2 — -  \  dr  r       e        \l-e     / 

1         2z+a   i=0  (V-DKZ+?)1  Jo   2  2 


x   sin(kr2    +  ~-  ^   2kr2+  ^ ) 
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Finally,  making  the  variable  change  x  =  (Z  +  d)r 


A1(V,X,2)  = 


32n2vl 


(2Z-KV)  (Z+a) 


^TlL^")1  (^ITT  Ii(v+x-i'  <w 


where 


i(y)=S ( 


x  il-e 


JXx  \  3 


(Z+a)\       .    /kx     Z-l  „  2kx     >,  .._. 

v    '  J  sml +  —r-  2n +  77/dx        (49) 


when  Z  =2  (He  ),  equation  (49)  must  be  evaluated  numerically.   How- 
ever when  Z  -  1 ,77  =0  and  the  integral  reduces  to  a  standard  form 
r3°l •   Numerical  evaluation  of  (49)  is  discussed  in  Appendix  D. 


Figure  2.   Interchanging  the  r   and  r   Integration 


Proceeding  now  to  the  evaluation  of  A„(v,X,2),  note  that 
the  r   integration  may  be  done  directly  and  that  evaluating  it  at 
the  ufjper  limit  (*-,)  gives  the  same  function  of  r   as  did  the 
evaluation  of  the  r   integration  at  the  lower  limit  (r  )  in  the 
evaluation  of  A, (v,X>2)  but  with  v  and  X  interchanged: 
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a      N  3      a. 


^     ^       32n2X  .'     (7    2      YV"         (\        ~2ri)       "2*1    V     .     ,,  Z-l0     „, 


A?(V 

~0 


v  m       ,  a      \3 


V      (2      \i         1  C°\      (\       "2ri)       "(Z+a')ri    V+A-i     .     ,.  Z-l„    01  „,1 

-    >       I- J      Tn — tt-t  \    drVl-e  /    e  r,  sin(kr,+ — ~-Sm  2kr    +77  )  \ 

L      \2Z+a/       (X-i).'    J         1  1  v       1         k  1     '  ') 

i=0 

Making    the    variable   change   x   =  Gtr    /2    in    the   first    term  and 
x   =    (Z+Gi)r,    in    the   remaining    terms,    and   recognizing    that    the 
summation    is    in    fact    A.(X,V,2)    gives 

A2(V,X,2)=      ***  A  '    ; J-   I     (V)-    A    (X,V,2)  (50) 

2  (2z-K*)av   1(2z+a)A     s 


where 

00  V 

T    /.  x      C      "x      y,-.       -xv3    .    f2kx         z-1  „     4kx  V  . -_  . 

Is(y)=  ^    e         x'(l-e      )    sin^-—   +  -—  ^  ^--   +  ry^dx    .  (51) 

Now  since  A,  (V,X,2)  and  A„(A,v,2)  always  occur  together  it  is  con- 
venient to  define 

AS(X,V,2)  =  A1(V,X,2)  +  A2(X,V,2)  (52) 

and  similarly 

As(X,v,u)  =A1(V,X»^)  +A2(X,v,M.)  (53a) 

Bs(X,v,^)  =  B1(v,X,n)  +  B2(X,v,n)  (53b) 

Ac(\,v5u)  =  A(v,\,[i.)    +A,(A,v,M.)  (53c) 

Bc(A,V,H)  =  B3(V,X,^)  +  BZt(X,v,M-)  (53d) 

where  (53a)  and  (53b)  arise  in  evaluating  <  s|  ( H-E) |x  >  and  (53c)  and 

(53d)  arise  in  evaluating  <c|(H-E)|X>   and  the  number  of  terms  in 

each  is  now  reduced  to  26.   Unfortunately,  explicit  expressions  for 

A  (v,X,2)  and  A„(v,X,2)  are  required  in  order  to  evaluate  A  (v,X,l), 
J-  o  s 
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B    (v,X,l),    A    (v,X,l)    and    B    (v,X,l),    respectively    so   the   reduction 

in    labor    is    not   so   great    as    it   might    appear.       Note    that   the   factor 

2 

32n    /(2Z+0!)    is    common    to    every    term,    so    it    may    be   dropped   from   the 

calculation . 

3.   Summary  of  the  Harris  Method 

To  summarize  the  Harris  prescription  as  it  applies  to 
atomic  scattering,  the  computational  scheme  is  as  follows: 

a.  For  the  given  basis  set  size  evaluate  the  required 
"close-in"  integrals  from 

A(vA;2)=x;{v:._l..IX?=o  4f!4^iLL}  ,54, 

b.  Evaluate  the  required  A(v,X,ia)  and  B(V,X>^)  using  the 
recursion  relations  of  Appendix  B. 

c.  For  a  given  CL ,  evaluate  H.  .  and  L.  .  using  the  expressions 

V  +X  +LL 

in  Appendix  A.   Note  that  each  A  and  B  must  be  divided  by  Of 
to  give  the  correct  result. 

d.  Solve  the  eigenvalue  problem  (24)  for  the  eigenvalues 
E  and  corresponding  eigenvectors  c  . 

e.  Find  an  appropriate  incident  electron  momentum  from 
the  expression 

kp  =  2Ep  +  z2  (55) 

2 
where  -Z   gives  the  ground  state  energy  of  the  scattering  center 

in  rydbergs  (1  ryd  =  13.6  eV) . 

f.  For  this  value  of  k  and  the  previously  assigned  value 

of  a,    evaluate  (numerically  if  necessary)  I, ,  I  ,  I   and  I   for  the 

'      1   s   3      c 

required  values  of  the  index  y   using  (49)  and  (51)  for  I   and  I   and 
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"S^ 


^     /     x      r       -*      >Y,  2(Z+0;)    V         /    kx         Z-l  „     2kx  \.  /c... 

x3(y)ssv    x A1_e  J  cosVi^ +  ir^i^ +T7;dx     <56> 

and 

CO 

Ic(y)=  3   e        Jf  (l-e      )    cos    l^—   +  —  &z  —   +  ^jdx  (57) 

g.       Using    these    values,    evaluate    the   required    terms    of   A    , 

A    ,    A      and   A„    for    u-    =    2    using    (1+8)    and    (50). 

h.       Using    the   recursion    relations    in   Appendix    B  evaluate 

the   remaining    required   terms    of   A    ,    B    ,    A      and    B    . 

s    s    c       c 

i.   Evaluate  the  <    s| (H-E) lx •  >  and  <  c| (H-E) |x •  >  using 

N 
the  expression  in  Appendix  A  and  form  £   c.  <S  |  ( H-E)  |X  •  >  ar>d 

N 

E   c  <S|(H-E)|X-  >  • 
i=l  l0  x 

j.       Finally  the  negative  quotient  of  the  two  sums  above  is 

tan6(k    ) ,    the    desired   result,    and   repetition    of   the  above    step    for 
r 

various   0!    normally   permits    the    entire    desired    range   of    k      to   be 

1  ^  p 

cover  ed. 
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III.   RESULTS  AND  DISCUSSION 

The  calculations  reported  in  this  section  were  performed  in 
double-precision  ( 16  figures)  arithmetic  on  the  IBM  360/67  computer 
of  the  Naval  Postgraduate  School  Computer  Facility.   The  programming 
language  was  Fortran  IV  and  the  H-level  compiler  was  used.   The  sub- 
routines used  to  factor  the  matrix  L  and  invert  U  and  to  solve  the 
eigenvalue  problem,  and  the  32-point  Gauss -Legendre  quadrature 
subroutine,  as  well  as  several  auxiliary  subroutines  were  furnished 
by  the  source  library  of  the  Computer  Facility. 

The  principal  program  was,  of  course,  written  to  evaluate  the 
phase  shifts.   This  program  and  its  required  subroutines  are  re- 
produced following  the  appendices.   In  addition  it  was  found  con- 
venient to  have  available  modifications  of  the  main  program  which 
solved  only  the  eigenvalue  problem,  giving  in  one  case  the  values 
of  the  incident  electron  momenta  in  the  elastic  scattering  range 
and  in  the  other  the  calculated  energy  levels  of  the  bound  state 
system.   Since  the  scheme  adopted  for  checking  the  eigenvalues  and 
eigenvectors  was  rather  time  consuming  and  involved  a  rather  large 
quantity  of  printed  output,  it  was  done  on  only  selected  matrices 
and  then  only  in  the  supplementary  programs.   The  modification  of 
the  phase  shift  program  to  perform  these  supplementary  functions  is 
obvious  and  the  programs  so  modified  are  not  reproduced  here. 
Several  features  were  built  into  the  programs  to  permit  the  maximum 
flexibility  with  the  minimum  of  internal  changes  and  as  it  finally 
evolved  only  the  storage  capacity  within  the  program  had  to  be 
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changed,  depending  on  the  size  of  the  basis  set  to  be  used.   All 
other  options  were  controlled  by  data  inputs  from  external  devices. 

The  results  for  e-H  scattering  will  be  presented  first,  followed 

+ 
in  part  B  by  the  e-He   results.   Comparison  with  the  results  of 

other  calculations  will  be  made  where  appropriate.   The  results 

will  be  presented  in  the  following  order:  (i)  phase  shifts  for 

singlet  and  triplet  scattering;  (ii)  location  of  resonance  levels; 

and  (iii)  bound  state  energy  levels  from  the  eigenvalue  problem 
(see  Appendix  C) . 

A.   ELECTRON -HYDROGEN  SCATTERING 
1 .   Phase  Shifts 

Singlet  and  triplet  phase  shifts  were  calculated  for  basis 
sets  from  N  =  22  elements  (M  =  4,  singlet;  M  =  5 ,  triplet)  to 
N  =  70  elements  (M  =  7,  singlet;  M  =  8,  triplet).   Although  the 
program  had  built  into  it  the  capability  of  calculating  phase 
shifts  for  the  95-element  basis  set  (M  =  8,  singlet)  the  eigen- 
value problem  began  to  show  signs  of  instability  at  this  size  and 
the  computer  time  required  to  find  eigenvalues  and  eigenvectors 
became  unacceptably  long  (8-10  minutes  per  matrix),  permitting  no 
more  than  testing  the  program  at  this  level.   In  addition  the 
storage  area  required  for  such  programs  approached  the  limit  avail- 
able on  the  computer . 

Preliminary  to  calculating  the  phase  shifts,  equation  (24) 
was  solved  for  a  wide  range  of  values  of  the  non-linear  parameter  (X 
and  the  resulting  values  of  incident  electron  wave  number  were 
plotted  as  a  function  of  a. .   As  expected,  the  eigenvalues  fell  into 
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clearly  discernable  trajectories  whose  evolution  as  N  was  increased 
was  easily  traceable,  at  least  in  those  regions  where  the  eigen- 
values were  sparse.   Figure  3  shows  such  a  plot  for  the  singlet 
states  of  hydrogen  with  N  =  70 . 

This  plot  clearly  illustrates  a  major  disadvantage  of  the 
Harris  method,  which  has  been  alluded  to  by  others  Ll3>25l  in  the 
context  of  the  inconvenience  of  being  unable  to  pick  the  scattering 
energy  without  the  expenditure  of  considerable  labor.   The  dis- 
advantage seems  to  be  of  a  more  fundamental  nature,  however.   While 
coverage  of  the  entire  energy  range  is  certainly  possible,  at  least 
in  most  cases  (it  was  not  possible  with  the  basis  sets  used  to  reach 
zero  energy  in  the  e-H  triplet  case),  sufficient  points  at  any  given 
value  of  k  are  not  always  available  to  permit  the  sort  of  investi- 
gation of  convergence  rate  that  characterized  the  analyses  of 
Schwartz  L 12 J  and  Armstead  L  31 J  -   Furthermore  in  those  regions 
where  few  eigenvalues  exist  there  is  no  guarantee  that  they  will 
exist  for  an  optimum  value  of  a. .      Schwartz  has  indicated  that  the 
useable  range  of  0L   was  between  about  0.8  and  4.0  in  his  calculation, 
and  this  generally  proved  true  in  the  present  case  also.   Reference 
to  Figure  3  shows  that  even  in  the  largest  basis  set  used  the  value 
of  a  for  the  only  trajectory  passing  k  =  0.1  is  around  5.0,  just 
outside  this  range.   In  this  region,  increasing  N  is  not  likely  to 
produce  additional  eigenvalues  at  which  to  solve  the  phase  shift 
problem,  for  as  is  shown  in  Appendix  C,  the  eigenvalues  result  from 
a  Ray leigh-Ritz  variational  calculation  of  the  energy  levels  on  the 
two-electron  bound  system  (see  part  3  of  this  section)  and  hence 
the  eigenvalue  trajectories  shown  in  Figure  3  represent  approximations 
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Figure  3.  Eigenvalue  Trajectories  for  Singlet  e-H  System  N  =  70 
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to  the  various  energy  levels  of  H   for  a  given  OL.      Although  the 
scale  is  not  the  energy  scale  one  usually  associates  with  such 
plots  the  relationship  between  the  two  is 

k2  =  2E  +  1  (58) 

in  the  case  of  hydrogen,  where  E  is  given  in  atomic  units  (27.2  eV) . 
In  such  calculations,  as  N  increases  the  lowest  lying  trajectories 
become  quite  good  approximations  to  the  energies  of  the  lowest  lying 
eigenstates  of  the  bound  system  [32]  after  which  further  increase 
in  N  will' give  rise  to  no  additional  trajectories  in  that  region. 
In  fact  the  situation  is  likely  to  get  worse  from  the  point  of  view 
of  the  scattering  calculation,  since  as  N  is  increased  the  minima 
in  the  trajectories  become  broader,  approaching  horizontal  straight 
lines  as  completeness  is  approached.   Consider  an  incident  electron 
wave  number  near  but  above  that  corresponding  to  the  highest  well 
defined  energy  level  of  the  bound  system  for  a  given  N.   As  N  is 
further  increased  the  values  of  OL   for  which  such  a  wave  number  is 
an  eigcnstate  of  the  system  will  be  pushed  closer  to,  and  perhaps 
beyond,  the  limits  of  the  range  in  OL    for  which  reliable  results  can 
be  expected.   Figure  4  illustrates  this  behavior  for  a  typical 
eigenvalue  trajectory.   As  N  is  increased  the  curve  becomes  broader 
and  its  minimum  approaches  a  limit.   For  N  >  70  the  minimum  probably 
will  not  decrease  appreciably  but  the  curve  will  become  broader. 

One  might  conclude  that  in  an  energy  region  which  is  rich 
in  eigenvalues  the  situation  ought  to  improve.   That  such  is  not 
always  the  case  can  be  seen  by  examining  Figures  5  and  6,  showing 
singlet  phase  shifts  for  k  =  0.4  and  0.8  plotted  as  a  function  of  ft. 
It  must  be  noted  that  the  lines  connecting  the  points  arising  from 
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Figure  4.  Typical  Evolution  of  One  Eigenvalue  Trajectory  as  N 
is  Increased. 
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the  calculations  for  a  given  basis  set  are  just  that  and  no  more. 
They  have  no  physical  or  mathematical  meaning, for  the  points 
plotted  represent  a  complete  collection  of  the  possible  values 
to  be  found  at  a  given  energy  (possibly  neglecting  one  or  two 
at  small  a)  . 

In  Figure  5,  while  the  "curves"  of  all  basis  sets  have 
clear  cut  minima,  there  is  no  clearly  defined  trend  as  the  basis 
set  size  increases,  as  Schwartz  and  Armstead  observed  in  their  Kohn 
calculations.   This  may  be  due  to  numerical  errors  which  accumulate 
more  rapidly  as  the  basis  set  is  increased  in  size,  however  as  will 
be  seen  below  the  ambiguity  is  considerably  reduced  in  the  triplet 
case  which  leads  one  to  conclude  that  this  error  is  probably  small. 
Note  however  that  the  range  of  Oi   over  which  the  phase  shift  is 
nearly  constant  is  greatest  for  the  case  N  =  70  and  hence  the 
minimum  of  that  curve  has  been  selected  as  the  most  probable  value 
in  this  calculation.   The  rather  large  error  indications  in  Table 
II  reflect  this  situation. 

Referring  now  to  Figure  6,  note  that  while  the  minimum  of 
tHe  N  =  70  curve  is  at  6  =  0.76l,  there  is  a  plateau  forming 
around  6  =  0.886.   For  this  value  of  k  increasing  N  may  well  bring 
improved  results  since  the  minima  of  more  trajectories  may  be  ex- 
pected to  move  into  the  region,  however  the  ambiguity  will  remain 
if  the  plateau  becomes  broader  and  stable  in  value.   A  similar 
structure  is  evident  in  the  calculations  for  k  =  0.7  and  0.6  as  well 
It  is  interesting  to  note  that  these  plateaus  correspond  to  the 
values  for  the  phase  shift  given  by  Schwartz  and  the  difference 
between  these  plateaus  and  the  minima  of  the  curves  increases  with 
increasing  k  until  at  k  =  0.866  the  difference  is  in  excess  of  10%. 
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Figure  5.  Singlet  Phase  Shift  for  e-H  as  a  Function  of  a 
k  =  0.4. 
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Although   they    do   not    discuss    it,    close    examination    of   the    results 
of    the   Harris    calculation    of   e-H   singlet   phase   shifts    by    Oberoi 
and   Callaway    L 1 3 J    using    the    close-coupling   approximation    seems    to 
indicate   that    their    results    may   be   as    much  as    4-5%   less    than 
Schwartz's    at    k    =   0.8,    although   their    graphical   presentation    of 
results    makes    such   a    comparison    difficult. 

Appeal    to   experiment    to   resolve   the   ambiguity    is    fruitless 
at    this    time    for    the   difference    in    total   scattering   cr oss -section 
at    k   =  0.8   for    the    two    values    of    6,    0.886  and   0.76l    is    only   about 
3%,    far    less    than    the   accuracy    of   any    experiments    conducted    to 
date   [33].      The   experiments    report    results    about    15%   less    than 
those   of   Schwartz   and  Armstead  at    the   upper    energy    limit    and   it    is 
worth    noting    that    the   present    result    is    also    lower    than    that    of 
Schwartz's.       However    because   of    the    large   experimental    errors    re- 
ported   (^    20%)    the   significance    of    this    fact,    if  any,    cannot   be 
estimated . 

Based   on    the    results    available    it    seems    more    wishful    than 
logical    to   ascribe    more    significance    to    the  plateaus    in    the   N    =    70 
curves    than    to   the   minima    and    in    the   absence    of   the   results    of 
Schwartz    the   lower    value  would   most    likely    have  been    chosen   as    the 
most   probable   value.       Therefore   the    minima   will    be   chosen    for 
cons  istency . 

The   extrema    of    the    curves    for    phase    shift    vs.    0L   generally 
are   minima    except    at    the   lowest    energies    calculated   and   they    seem 
to   trend  downward,    so    in    most    cases    the   value    of   the  phase    shift 
for    N    =    70   can   probably    be   taken   as    an    upper    limit    except    for 
k  ^    0.2    where    it    seems    to   be   a    lower    limit.       Where   an    unambiguous 
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Figure  6.  Singlet  Phase  Shift  for  e-H  as  a  Function  of  d 
k  =  0.8. 
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convergence  could  be  established  it  seemed  to  be  about  as  the 
sequence  1/N. 

The  scattering  length  is  defined  to  be 

tan5 
a   =  lira  — r— —     .  (59) 

k-0 

In  the  present  calculation  it  seemed  natural  to  evaluate  the 
scattering  length  for  e-H  scattering  by  simply  evaluating  (tan5)/k 
for  a  succession  of  small  values  of  k  and  extrapolating  the  results 
to  k  =  0.   This  was  done  for  the  case  N  =  70  and  the  smallest  value 
of  k  used  was  .000412.   The  extrapolation  to  k  =  0  gave 

a   =  -6.02181(1) 
o  v  ' 

where  the  figure  in  parentheses  is  the  uncertainty  in  the  last 
digit.   This  value  differs  by  about  1%   from  the  value  a   =  -5-965 
reported  by  Schwartz  [l2j.   Note  that  the  uncertainty  reported 
above  is  the  uncertainty  arising  from  the  extrapolation  to  zero 
and  is  not  necessarily  an  indication  of  absolute  accuracy.   The 
extrapolation  used  in  evaluating  a   is  shown  in  Fiaure  7. 
It  should  be  noted  that  the  difficulties  with  the 
Harris  method  discussed  above  should  not  be  peculiar  to  this  choice 
of  trial  wave  function,  for  since  the  Hamiltonian  (30)  is  the  exact 
Hamiltonian  for  any  two-electron  system  the  eigenvalue  calculation 
should  lead  to  the  same  effects  described  above  for  any  choice  of 
trial  wave  function  as  long  as  it  meets  the  usual  boundary  condi- 
tions and  is  based  on  a  quadratically  integrable  basis.   Perhaps 
by  using  some  approximate  Hamiltonian  the  eigenvalue  situation 
would  be  improved  but  the  loss  of  accuracy  from  this  approximation 
would  probably  offset  any  gains  from  improved  eigenvalue  behavior. 
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Figure  7.  Evaluation  of  e-H  Singlet  Scattering  Length 
N=70 
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The  results  for  the  triplet  phase  shifts  are  much  as  has 
been  described  above  for  the  singlet  case  except  that  the  curves 
seem  to  trend  upward  toward  maxima  for  values  of  k  larger  than  was 
found  in  the  singlet  case.   Since  in  the  triplet  case  the  electron 
spins  are  parallel  the  effective  repulsion  between  the  electrons 
due  to  the  Pauli  principle  tends  to  improve  the  convergence.   Thus, 
while  the  qualitative  picture  is  similar  to  that  of  the  singlet 
scattering,  the  ambiguities  are  much  less  severe.   Figure  8  shows 
the  eigenvalue  trajectories  for  N  =  70.   Note  that  there  is  no  tra- 
jectory that  reaches  k  =  0,  the  lowest  one  just  barely  passing 
k  =  0.1.   Hence  it  was  not  possible  to  evaluate  the  triplet  scat- 
tering length  as  was  done  for  the  singlet  case.   Figures  9  and  10 
show  plots  of  phase  shift  vs.  0i   for  k  =  0.4  and  0.8.   Finally,  the 
curve  of  6  vs .  k  for  both  singlet  and  triplet  scattering  is  shown 
in  Figure  11.   The  resonances  shown  on  the  plot  are  discussed  in 
the  next  section.   In  each  case  one  resonance  which  has  been  re- 
ported elsewhere  as  being  just  below  the  excitation  threshold  was 
observed  in  this  calculation  to  be  just  above  and  is  so  indicated 
on  the  plot  by  the  vertical  dashed  lines.   It  is  believed  that  in- 
creasing the  basis  set  size  would  bring  these  values  below  the 
threshold.   The  singlet  and  triplet  phase  shifts  are  tabulated  and 
compared  with  those  of  Schwartz  [l2]  in  Table  II. 
2 .   Scattering  Resonances 

The  present  calculation  offers  a  feature  not  found  in 
Schwartz's  calculat ion--the  localization  of  real  scattering  re- 
sonances.  There  is  probably  no  theoretical  reason  why  the  Kohn 
method  should  not  show  the  resonances,  however  since  they  are  very 
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Figure  9.  Triplet  Phase  Shift  for  e-H  as  a  Function  of  a 
k  =  0.4 
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Figure    11.    Calculated  Phase   Shifts    for    e-H   Scatter in< 
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narrow   [l8,34,35]    it    is    quite    likely    that    they   simply   were    not 
resolved.       It    might    also   be    difficult    to    distinguish    them    from    the 
spurious    ones    characteristic    of   such    calculations. 


Incident 
Electr  on 
Wave   Number 


Scattering    Phase    Shift    in    Radians 


0.1 
0.2 

0.3 
0.4 
0.5 

0.6 

0.7 
0.8 

0.9 


Singlet 


Present 
Work 

2.566(5 
2.067(4 
1.696(6 

i.4i4(6 
1.194(7 
1.020(4 

.878(3 
.763(4 
.693(7 


Schwartz[l2] 

2.553 

2.0673 

I.6964 

1.4146 

1.202 
1.04l 

.930 

.886 


Tr  iplet 


Present 
Work 

2.9386(5) 
2.7187(7) 
2.5013(7) 
2.2980(5) 
2.1093(3) 
1.9321(5) 
1.7772(4) 
1.6389(5) 
1.5622(6) 


Schwartz[l2] 

2.9388 

2.7171 

2.4496 

2.2938 

2 . 1046 

1.9329 

1.7797 

1.643 

1.563 


The    figures    in   parentheses    represent    the   uncertainty    in    the   last 
digit    reported. 

Table    II.       Singlet    and  Triplet   Phase    Shifts    for    e-H 
Elastic    Scattering. 


The    existence    of    scattering    resonances    has    long   been 
associated   with    the  auto-ionizing    states    of    the    compound   system   [36] 
and   recent    experiments    and   calculations    have    confirmed   this 
L35s37,38].       Auto-ionizing    states    (in    the    case    of   H      they   are   more 
properly    called   auto-detaching)    are    those    in   which    the    total    energy 
of    the    system   exceeds    the    ground   state    energy    by   more    than    the 
ionization    energy    of   the    two-electron    system    but    in    which   both 
electrons    are    in    excited   states.       Such   states    can    decay    non- 
radiatively    when    one   electron's    excess    energy    is    transferred    to 
the    other    which    is    then    ejected   from    the  atom   while   the   first 
reverts    to   the    one-electron    ground   state.       Such    states    are   normally 
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extremely  short  lived,  typical  lifetimes  often  being  several  orders 
of  magnitude  less  than  the  normally  allowed  radiative  transitions  L39l 
Clearly  an  electron  captured  momentarily  into  such  a  state  will  be 
considerably  delayed  in  its  passage  by  the  atom,  suffering  a  cor- 
respondingly greater  phase  shift,  and  thus  show  a  resonant  behavior 
in  the  cross  section.   Since  these  states  lie  above  the  one  electron 
ionization  potential  and  are  extremely  short-lived  they  are  often 
referred  to  as  quasi-bound  [4o]. 

Figures  3  and  8  show  the  evidence  of  such  states  in  the 
singlet  and  triplet  structure  of  H  .   They  are  the  broad  shallowly 
curved  trajectories  lying  just  below  the  excitation  threshold  in 
each  case. 

Since  these  states  do  lie  above  the  ionization  potential 
of  H   and  thus  among  a  continuum  of  states,  one  expects  from  the 
Hylleraas-Undheim  theorem  [32]  that  such  states  would  not  be  well 
approximated  by  an  eigenvalue  trajectory  until  after  all  lower 
lying  states  had  been.   The  reasons  for  this  rot  being  so  are  not 
well  understood,  however  Hol^ien  and  Midtal  C40]  have  found  that  if 
the  trial  wave  functions  include  an  appropriate  mix  of  the  one 
electron  state  functions  (i.e.,  lsns  ,  n  >  1,  and  also  Isls,  2s2s  , 
2p2p,  etc.)  that  the  eigenvalues  corresponding  to  the  autoionizing 
states  begin  stabilizing  long  before  those  of  many  lower  lying 
states.   They  conjecture  that  this  probably  occurs  because  the 
equal  treatment  of  the  two  electrons  leads  to  an  approximate 
or thogonal ization  to  the  lower  states.   The  limits  approached  by 
these  trajectories  may  not  be  upper  bounds  however. 
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That  such  conditions  occur  in  the  present  problem  is  most 
likely  a  result  of  the  use  of  Hylleraas-type  trial  functions,  which 
can  be  combined  to  approximate  the  one-electron  orbitals  needed  to 
satisfy  the  condition  set  by  HoljtSien  and  Midtal  above. 

In  the  present  calculation  two  resonances  were  resolved 
within  the  e-H  elastic  scattering  range.   One  in  the  singlet  states 
was  observed  at  k  =  O.83713  which  corresponds  to  a  bound  state 
energy  of  E  =  -.149610  atomic  units  (1  a.u.  =  27.2  eV) .   A  second 
in  the  triplet  states  was  observed  at  k  =  0.864l8  which  corresponds 
to  E  -  -.126596  a.u.   In  addition  there  were  two  such  states  identi- 
fied just  above  the  excitation  threshold  which  probably  correspond 
to  levels  predicted  just  below  the  threshold  by  others  L 3§ J •   In 
all  liklihood  a  larger  basis  than  those  used  here  would  give  a 
better  approach  to  the  results  reported  by  others.   Since  the  means 
of  locating  these  resonances  was  from  a  bound  state  calculation  no 
estimate  of  their  width  can  be  given.   A  characteristic  of  these 
calculations  is  that  the  phase  shifts  associated  with  a  particular 
eigenvalue  trajectory  tend  to  be  rather  independent.   Thus  evaluating 
a  phase  shift  at  a  value  of  k  corresponding  to  a  resonance  energy 
but  at  a  value  of  GL   corresponding  to  a  normal  (i.e.,  non-resonant) 
eigenvalue  trajectory  appears  to  show  little  or  no  evidence  of  the 
resonance.   Although  the  non-resonant  and  resonant  trajectories 
sometimes  do  cross,  no  conscious  attempt  to  evaluate  the  phase 
shift  at  such  a  crossing  point  was  made  since  in  the  eigenvalue 
method  used  (Jacobi  variable  threshold  C 4l 1  the  eigenvectors  found 
for  degenerate  or  nearly  degenerate  eigenvalues  may  be  grossly  in- 
accurate (Ref.  4l,  p.  28l)  and  thus  the  phase  shift  calculated  at 
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this  point  would  be  difficult  to  interpret.   The  independence  of 
the  phase  shifts  corresponding  to  different  trajectories  is  probably 
related  to  the  orthogonality  properties  inherent  in  the  eigenvalue 
problem  and  discussed  above  but  beyond  this  it  is  not  well  understood 
Table  III  summarizes  the  results  for  the  e-H  resonances  and  compares 
them  with  other  calculations  and  with  such  experimental  data  as  is 
available. 

3.   H~  Energy  Levels 

As  shown  in  Appendix  C  and  as  demonstrated  in  the  last 
section,  the  first  part  of  the  Harris  prescription  involves  what 
amounts  to  a  variational  calculation  of  the  bound  state  energies. 
Reports  of  such  calculations  abound  in  the  literature,  commencing 
ith  Hylleraas  C20]  and  continuing  to  the  present,  and  have 
ttai'ned  a  high  degree  of  accuracy.   Hence  the  results  obtained 
here  are  of  no  more  than  passing  interest  except  that  they  arise  as 
a  natural  by-product  of  the  Harris  method,  and  this  fact  seems  not 
to  have  been  mentioned  previously  in  the  literature  concerning  the 
Harris  method.   Considering  that  no  special  pains  were  taken  to 
assure  accuracy  (other  than  the  normal  and  reasonable  ones  dis- 
cussed in  Appendix  D)  the  energy  results  are  surprisingly  accurate, 
and  since  they  are  obtained  with  no  additional  effort  from  the 
phase  shift  calculation  it  seems  worthwhile  to  discuss  them  briefly 
here. 

Peker is  C46]  has  reported  a  calculation  of  the  ground  state 
energies  of  He  and  H   in  which  he  used  matrices  of  up  to  order  IO78 
and  was  at  great  pains  to  achieve  the  maximum  possible  accuracy. 
He  reports  the  ground  state  of  H   from  this  calculation  to  be 
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-.52775097  a.u.   In  the  present  work  a  value  of  -.5277477  a.u.  was 
obtained  with  an  order  70  matrix. 

As  an  example  of  the  bound  state  results  obtainable, 
Figures  12  and  13  show  the  negative  energy  structure  for  the 
singlet  (N  =  70)  and  triplet  (N  =  50)  calculations.   Not  all  the 
detail  is  included  in  the  region  near  zero  energy  because  the  large 
number  of  eigenvalues  in  that  region  makes  it  very  difficult  to 
sort  out  the  eigenvalue  trajectories.   Note  that  the  lower  eigen- 
values clearly  exhibit  the  broad  minima  alluded  to  in  section  1. 
This  characteristic  of  the  trajectories  is  even  more  pronounced  in 
the  case  of  helium,  which  will  now  be  discussed. 

B.   ELECTRON -HELIUM  ION  SCATTERING 
1 .   Phase  Shifts 

In  principle  the  only  difference  between  scattering  elec- 
trons from  helium  ions  and  from  hydrogen  is  the  presence  of  the 
long  range  Coulomb  force  in  the  former.   This  effect  is  accounted 
for  in  the  form  chosen  for  the  asymptotic  part  of  the  trial  wave 
function  and  once  the  more  difficult  problem  of  numerically  evalu- 
ating the  resultant  integrals  is  solved,  one  expects  the  procedure 
to  be  straightforward  and  to  yield  results  whose  accuracy  is  de- 
graded only  by  the  relative  inaccuracy  of  numerical  integration  vs. 
exact  integration.   Appendix  D  contains  extensive  discussion  of  the 
method  of  evaluating  the  numerical  integrals  and  assuring  their 
accuracy,  and  it  appears  that  a  high  degree  of  confidence  can  be 
placed  in  them.   Nevertheless  the  results  of  these  phase  shift 
calculations  are  not  so  accurate  as  the  hydrogen  ones,  if  one 
accepts  the  results  of  Burke  and  Taylor  Cl8]  as  the  most  accurate 
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Figure  12.  Singlet  Energy  Levels  for  H~.  N  =  70 
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Figure  13.  Triplet  Energy  Levels  for  H  .  N  =  50 
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to  date.   The  close  agreement  of  their  results  with  those  of 
others  by  different  methods  [35,44]  indicates  that  this  confidence 
is  probably  well  placed.   Thus  it  would  appear  that  the  Coulomb 
force  has  an  effect  on  the  accuracy  of  the  calculation  beyond  what 
is  initially  expected. 

As  Figures  14  and  17  show,  the  eigenvalue  trajectory  situ- 
ation is  somewhat  improved  in  the  case  of  helium,  owing  to  the 
existence  of  a  complete  set  of  two  electron  bound  states  in  the 
region  below  the  first  ionization  potential  of  He.   Nevertheless 
the  singlet  phase  shifts  at  low  energy  appear  to  be  as  much  as  an 
order  of  magnitude  larger  than  the  results  reported  by  Burke  and 
Taylor .   Because  the  calculation  gives  only  tan6  the  value  of  5 
is  uncertain  by  +_  n  ^ ,  and  the  choice  of  6  to  be  .^  at  k  =  0.2  is  in 
closer  agreement  with  Burke  and  Taylor.   However  the  slope  of  the 
6  vs .  k  curve  at  k  =  0.2  was  measured  and  its  value  seems  to  support 
the  higher  value  chosen.   At  higher  energy  the  convergence  appears 
to  be  to  a  result  about  5-10rc  higher  than  that  reported  by  Burke 
and  Taylor,  although  the  existence  of  the  scattering  resonances 
at  higher  energies  undoubtedly  affects  the  results  and  makes  com- 
parison difficult  in  this  region.   The  existence  of  the  resonances 
makes  more  likely  the  possibility  of  attempting  to  evaluate  the 
phase  shift  at  a  nearly  degenerate  eigenvalue,  with  its  concomitant 
inaccurate  eigenvector  which  may  further  confuse  the  picture. 
Typical  plots  of  5  vs .  0L   for  singlet  e-He   scattering  are  shown  in 
Figures  15  and  16 . 

A  similar  but  (not  unexpectedly)  less  severe  situation 
exists  with  the  triplet  scattering.   Figures  18  and  19  show  typical 
triplet  plots  of  6  vs.  0L . 
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Figure   l4 .    Eigenvalue  Trajectories    for    Singlet    e-He 
System.    N=70 
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Figure  15-  Singlet  Phase  Shift  for  e-He   as  a  Function  of  a. 
k  =  0.8 
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Figure  16 .  Singlet  Phase  Shift  for  e-He   as  a  Function  of  C£ 
k  =  1.4 
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Figure    17.    Eigenvalue   Trajectories    for    Triplet    e-He 
System.    N=50 
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Figure  18.  Triplet  Phase  Shift  for  e-He   as  a  Function  of  a 
k  =  0.8 
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a 

Figure  19.  Triplet  Phase  Shift  for  e-He   as  a  Function  of  a. 
k  =  1  .4 
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Results    of    the  present    calculation  are  presented  and   com- 
pared  with    those    of    Burke   and  Taylor    in   Table    IV   and  plotted   in 
Figure    20. 

2 .       Scattering   Resonances 

It   was   possible    to    identify   six    singlet    resonances   and   four 
triplet    resonances.       As    shown    in   Table   V    the   agreement   with   the   re- 
sults   of    Burke   and   McVicar    [44]    and   others    [37,40,42,43,471    is 
quite    good. 

Burke   and  McVicar    give   the    energies    of   seven   and   five   re- 
sonances   for    singlet   and   triplet    states,    respectively.      As    in    the 
hydrogen   case   an    additional    resonance    level    could   be    identified   just 
above    the    excitation    threshold  by    examining    the   energy    level   plots 
for    each    system.      Again    as    in    the   hydrogen,    these    levels    would 
probably   appear    just    below   the   excitation    threshold   if   the    calcu- 
lation   were    done    with   a    sufficiently    large   basis    set. 
3 •       Helium    Energy    Levels 

The    helium    energy    spectrum    has    been    even    more    extensively 
calculated   than    has    hydrogen    [20,46,48,49.].      The   most    accurate   such 
calculation    is    that    of  Pekeris    [46]    who  has    reported   a    value    of 
-2.903724375   a.u.    for    the   ground   state    energy    of   He.      After    applying 

mass   polarization    and   relativistic    corrections    and   correcting    for 

5         -1 
the   Lamb   shift,    he    finds   a    value    of   I.98310687   x    10      cm        for    the 

ground   state    of   He,    which    compares    with    the   experimental    deter- 
mination   by    Herzberg    [50]    of    I.983IO8    x    10      +_    .05    cm"    .       Similarly, 
for    the   2JS   state    of  He,    Pekeris    obtained   -2.17522937822   a.u.    or 
3.84566   x   lO^cm"      while   Herzberg   measured   3.845473   x    10^   +    .05   cm"    . 
In    the  present    calculation    -2.90372410    a.u.    and    -2.175191   a.u.    were 
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Scattering  Phase  Shift  in  Radians 


Incident 

Singlet 

Triplet 

Electron 

Present 

Burke    and 

Present 

Burke   and 

Wave   Number 

Work 

Taylor    [l8] 

U'or  k 

Taylor    [l8 J 

0.2 

3-5(9)a 

A- 1(9) 

0.4 

.98(4) 

1.469(8) 

0.4472 

.4217 

.9088 

0.6 

.63(2) 

1.108(6) 

0.6325 

.4092 

.8873 

0.7746 

.3989 

.8649 

0.8 

.512(8) 

.965(6) 

0.8944 

.3912 

.8462 

1.0 

.446(7) 

.3850 

.889(4) 

.8253 

1.1832 

.3780 

.7895 

1.2 

.417(7) 

.831(5) 

1.3416 

.3779 

.7593 

1.4 

.407(4) 

.778(4) 

1.4832 

.3910 

.7341 

1.5811 

.2977 

1.6 

.348(4) 

.730(4) 

1.6553 

.3988 

I.6712 

.7006 

1.6971 

.6931 

1.7073 

.2572 

1.732 

.326(3) 

.705(4) 

The  figures  in  parentheses  represent  the  uncertainty  in  the  last 
digit  rep or  ted. 


Table  IV.   Singlet  and  Triplet  Phase  Shifts  for  e-He 
Elastic  Scattering 
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Figure  20.  Calculated  Phase  Shifts  for  e-He   Scattering 
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1-3 

obtained  for  the  1  S  and  2  S  levels  in  He,  respectively.   Figures 

21  and  22  show  the  calculated  energy  levels  for  the  singlet  and 
triplet  S  states  of  He. 
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Figure  21.  Singlet  Energy  Levels  of  He.  N=70 
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Figure  22.  Triplet  Energy  Levels  of  He.  N=50 
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IV.   CONCLUSION 

The  principle  advantage  in  calculating  phase  shifts  by  the 
Harris  method  is  that  it  is  not  required  to  evaluate  matrix 
elements  involving  only  the  asymptotic  part  of  the  wave  function. 
This  results  in  a  considerable  computational  simplification  but 
the  price  paid  is  that  the  error  is  now  of  first  order  and  hence 
the  method  is  inherently  less  accurate  than  methods  dependent  on 
the  Kato  identity. 

An  additional  disadvantage  is  that  in  certain  regions  there  may 
be  a  paucity  of  energy  eigenvalues  at  which  to  carry  out  the  calcu- 
lation. Since  these  eigenvalues  are  approximations  to  the  compound 
system  energy  levels,  in  those  regions  where  the  approximations  are 
good,  increasing  the  number  of  parameters  in  the  trial  wave 
function  may  simply  make  the  situation  worse.  Even  in  regions  rich 
in  eigenvalues  the  fact  that  they  are  related  to  the  bound  states 
may  have  unpredictable  effects  on  the  phase  shift  calculations. 

The  Harris  method  does  have  a  feature  that  could  prove  useful 
in  complex  systems.   Resonances  in  the  scattering  cr oss -section 
can  be  uniquely  identified  as  arising  from  certain  bound  states  of 
the  compound  system,  such  as  the  auto-ionizing  states  in  the  pre- 
sent instance.   Since  the  scattering  at  a  particular  energy  can  be 
identified  with  the  compound  state  energy  level,  when  that  level  is 
well  defined  by  the  trial  wave  function  the  correlation  with  the 
behavior  of  the  phase  shifts  in  that  region  will  be  clear. 
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Whether  the  anomalous  results  for  the  e-He   system  found  here 
are  an  indication  of  a  defect  in  the  Harris  method  or  due  to  a 
subtle  error  in  computation  or  analysis  on  the  part  of  the  author 
is  not  absolutely  clear.   Although  the  latter  cannot  be  ruled  out, 
the  success  of  the  method  in  calculating  the  e-H  scattering  using 
the  same  computer  program  and  even,  as  pointed  out  in  Appendix  D, 
the  same  integration  scheme  (on  a  test  basis),  seems  strong  evi- 
dence in  favor  of  the  former  conclusions.   A  third  possibility,  of 
course,  is  that  the  results  obtained  here  are  correct  and  the  re- 
sults of  Burke  and  Taylor  and  others  are  incorrect.   In  the 
absence  of  experimental  evidence  to  the  contrary,  the  consistency 
of  the  other  calculations  would  seem  to  make  this  conclusion  also 
unlikely . 

Thus  it  is  the  tentative  conclusion  of  this  thesis,  that  the 
Harris  method  should  be  applied  to  electron-ion  scattering  with 
some  degree  of  caution. 
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APPENDIX  A 


Formulae   for    the   Various    Matrix   Elements 


1  .       H.  .    and   S  .  . 
_2J U 

The  matrix   elements  <   Y.IhIy    >  is    made   up    of   52    terms,    each    of 

1  j 


which    includes    an    integral    of   the    form 


A(V,X  ,u)  =   \dr    dr 


-  a(r    +r  _)         _   v     „    ,,    _ 
^     _-•  v     1       2'    V-2   A -2    u-2 


r         r         r 
1  2  12 


(A-l) 


or 


a/       >        I       xV     T      '   a(ri  +  r2}    V-2   X-2    u-2 
B(v,A  ,a)=   <^dr1dr2e  rx      r2      r 


12   cosG 


12  ' 


(A-2) 


The  trial  wave  function  is 


X  =  e 


a. 

2[I1+I2}    f    n  £  i  n\  m 

lrir2  irir2>i; 


(A-3) 


and  can  be  characterized  by  the  three  explonents  n,£,m.   Let  \. 

character  ized  by  n.,  £ . .  m.  and  V  .  by  n.,  I ..    m..   With  this 

ill  J  D         D         J 

notation    H. .    can    be   written 
1J 


be 


H..    =  -    (a/4)-A(n.+n.+2,4.+4.+2,m.+m.+2 

ij  /vij'ij'ij 

+    [^(n  .+    1)-Z    ]  -A(n.+    n  .+    l.£.+   £  .+    2,m.+    m.+2 

J  i  J  i         J  i         D 

+    [^(jG  .+    1)-Z    ]  -A(n.+    n  .+    2,£.+   i  .+    l,m.+   m.+2 

+  JgQ!m.'A(n.+   n  .+    2  ,£  .+   £  .+    3,m.+   m. 

+  ^m.'A(n.+    n  .+    3,^.+    ^  .+    2,m.+   m. 

J  i         j  l         j  l         j 

-    m.(m.+   £  .+    n  .+    1)-A(n.+    n  .+    2,£.+   £  .+    2,m.+    m. 
^n.(n.+    1)-A(n.+   n.,£.+   £  .+   2,m.+    m  .+   2 
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%£  .(4  .+   1) -A 
3       3 


V*ra  . 

3 

B 

%Xm 

■B 

ra  .n  . 
3    3 

'B 

m  .4  . 
3    3 

■B 

*2 

•A 

+    [^(4  .+    1)-Zj  -A 

3 

+  [JgO£(n  .+  i)-z]  -A 

+  ^ttm.-A 

+  %2m  .  •  A 

3 

-    m  .(m  .+    4  .+    n  .+    1 )  -A 
J      J  J  J 

hi  .u  .+  i)  -a 
J    J 

\n  .  (n  .+    1)  -A 


-'2»ra  . 
3 

•B 

B 

ra  .4  . 

B 

m  .n  . 
3    3 

'  B 

+    (V*2'A 


+    [%X(l  .+    l)-z]  "A 

+    [%a(n .+    l)-z]  -A 


n.+    n  .+  2,4.+   4.,m.+    m.+    2 
i         3  i         J       i         3 

n.+    n  .  +  2,4.+   4  .  +    2  ,  m  .  +   m.+    1 
13  13  13 

n.+    n.+  2,4.+    I  .+    3,m.+    m  . 
13  1313 

n.+    n  .+  3,4.+    £  .+    2,m.+   m. 
1313  13 

n.+   n  .+  1,4.+   4  .+    3,m.+   ra  . 
13  1313 

n.+   n  .+  3,4.+   4.+    l,m.+   ra. 
1313  13 

4.+   4  .  +  2,n.+   n  .  +    2  ,  m .  +   m.+    2 
13  13  13 

4.+   4.+  l,n.+    n  .+    2,m.+    ra.+    2 
13  1         3  13 

£.+    4.+  2,n.+    n  .+    l,m.+    m.+    2 
13  13  13 

£  .+   4  .+  2  ,n  .+    n  .+    3  ,ra  .+    ra  . 
13  1         j  1         j 

4.+   4.+  3,n-+   n  .+   2,m.+   ra. 
1        J  1        3  13 

4.+   t  .+  2  ,n  .  +    n.+    2  ,m.  +    ra. 
13  13  13 


4.+   4  .  ,n  .  +    n  .+    2,m.+    m.+   2 
1         J       1         D  13 

4.+    4.+    2,n.+    n.,ra.+    m  .+    2 
1  j  1  j       1         j 

4  .  +   4  .  +    2  ,  n  .  +   n.+    2  ,  ra  .  +   m.+    1 
13  1         3  13 

4.+    £  .+    2,n.+    n  .+    3,m.+   ra. 
13  1313 

4.+   4  .+   3,n.+    n  .+   2,ra.+   ra. 
1         3         '1         j         '1         j 

4.+   4.+    l,n.+    n  .+    3,m.+    m. 
x         j  1         j  1         3 

4.+  4.+    3  ,r> .  +   n.+    l,m.+   ra. 
13  13  13 

n.+   £  .+   2,4.+   n  .+    2,m.+   m .+    2 
13  13  13 

n.+   4.+    1,4.+    n.+    2,m.+    ra.+    2 
13  13  13 

n.+   4.+    2,4.+    n.+    l,m.+    ra.+    2 
13  13  13 
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+  J^ara  .  •  A 

J 

J 

-   m  .(ra  .+    n  .+   £  .+    1)  -A 
J      J         J         J 

%£  .(i  .+    1)  -A 
hn  .(n  .+    1)  -A 


^m  .  •  B 
^m  .  •  B 


ra  .£  . -B 
3    3 

ra  .n  .  "  B 
3    3 


+    [%3£(n  .+    1)-Z]  -A 

+  [^a(A  .+  i)-z]  -a 

J 

+  Jgttm  .  'A 

J 

+  ^am  .  *A 

3 

-    m  .  (ra  .+    n  .+    £  .+    1)  "A 
J       J         J         3 

^n  .(  n  .+    1)  -A 

J       3 

i>£  .(X  .+    1)  -A 
3      3 


■  A 


%2ra  .  •  B 
^m  .  -B 


m  .n  .  •  B 
3    J 


m  J  .  •  B 
3    J 


n.+   £  .+   2,£.+   n  .+    3,m.+   m. 
i         J  x        J  x        j 

n .  +   £  .+    3,£-  +    n  .+    2,m.+    m. 
l         j  l         j  l         j 

n .  +   £  .+    2,£.+   n  .+   2,m.+   ra. 

n.+   £.,£.+   n.+    2  ,ra.+  m.+    2 
l         3      l         3  x        j 

n  .  +   £.+2,£.+    n.,m.+   m  .+    2 
i         3  i         3      i         j 

n.+    £  .  +    2  ,  £  .  +    n.+    2  ,  m .  +   m.+    1 
13  13  13 

n . +    £  .+   2,£.+    n  .+    3,m.+   ra  . 
13  13  13 

n  .  +    X  .+    3J  .+    n  .+    2  ,ra.+   ra  . 
13  x         J  x         J 

n.+   £  .+    1,£.+    n  .+    3,m.+   m. 
^3  x        j  x        j 

n.+   £  .+    3  ,£  .  +    n.+    l,m.+   m. 
13  13  13 

£ .  +    n.+    2,n.+   £  .  +    2,m.+    m.+    2 
3  j  1         j  13 

£.+    n  .+    l,n.+   £  .+    2,ra.+   m  .+    2 
1         j  13  13 


£  .+    n.+  2,n.+   £  .+  l,m.+   ra  .+    2 
13              13  13 

£  .  +    n.+  2,n.+    £  .  +  3,m.+   m. 
13              1         3  J 3 

£.+    n  .+  3,n.+   £  .+  2,m.+   m  . 
13              x        j  '13 

£.+    n  .+  2,n.+    £  .+  2,m.+    m. 
13               13  13 

£.+   n.,n.+    £  .+   2,m.+    m  .+   2 

x        J  x        j  1        j 


.+    n  .+    2,n.+    £.,m.+    m.+    2 
13  x        3       1        j 


£.+    n  .+  2,n.+    £  .+    2,m.+   ra  .+    1 
13  13  1         3 

£.+    n  .+  2,n.+    £  .+    3,m.+   m. 
13  x        j  13 

£.+   n.+  3,n.+   £.+    2,m.+   m. 
x        J  x        j  13 

£.+    n  .+  l,n.+   £  .+    3,m.+   m. 
1         3  13  13 

£.+   n.+  3,n.+   £.+    l,m.+.  ra. 

1         3  •"    1         3         '1         3 


}  (a-;) 
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and 

S..=  A(n.+    n.+    2,i.+   £.+    2,ra.+   m.+    2) 

ij  1         j  1         j  1         j 


A(&.+   I  .+    2,n.+  n  .+    2,m.+   ra.+    2) 

+    {    A(n.+   £  .+   2,£.+  n  .+    2,m.+    in  .+   2) 

—             v     1          j               x  j               1          j          ' 

A(£ . +   n  .+   2,n.+  X  .+    2,m.+   m.+   2)    }     (A-5) 

v    l         j              l  j              l        j          '          v          ' 


Note    that    if   the    expression    (54)    is    used   to   evaluate    the  A's 

and    B's    that    each    of    the   terras    in    (A-4)    and    (A-5)    above   must    be 

v  +X  +J- 
divided   by  Oi 

2  *       <Sl  (H    -    E)  l\>    and  <C|  (H    -    E)|\  > 

Each    of   these   matrix   elements    is    made   up    of   26    terms    which    in 

turn    are   made   up    of   an    integral    of   the   form    (53    a-d)    where   A's    and 

B's    are   as    defined    in    (    44)    to    (    47)    and    immediately    following. 

Now    letting    \    in    (A-3)    above   be   characterized  by    n,£,m,   <S ,C |  (H-E)  \\  > 

is 

<S,C| (H-E)|X  >  =    ~  [H*   +    E) *A  (£+    1 ,n+    2,ra+   2) 

s  ,  c 

+    [^(4+    1)-ZJ-A  (4,n  +   2,m+    2  ) 

+   [^(n+    l)-z]-A  (1+    l,n  +    l,m+  2) 

+                           ^Q!m'A  {1+    l,n+    3, m  ) 

+                           ^m-A  (£+    2,n+    2,m  ) 

-ra(ra+   1+    n+    1)-A  {1+    l,n+    2,m  ) 

%£(£+    1)-A  (X-    l,n+    2,m+  2) 

^n(n  +    1)-A  (£+    l,n,m+    2  ) 

+                                     A  {1+    l,n+   2 ,m+  1) 
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-sara-B         (4  +    l,n  +    3,m 


^m-B         (4  +    2,n  +   2,m 


4m« B        (4 ,n+    3,m 
s  ,cv 

ran  *  B         (4  +   2  ,n  +    1  ,ra 
s  ,cv 


+    [-(%»    +    E)  -A         (n+    1,4+   2,m  +    2 


s  ,c 


+   [%0t(n+    l)-z]-A  (n,4  +   2,m+    2 


s  ,c 


+   [-ja(4+    l)-z]-A  (n+    1,4+    l,m+    2 

+  ^am«A  (n+    1,4+    3,m 

+  V*™  *  A  (n+2,4+2,m 

-ra(m+   4+   n  +    1)*A  (n+    1,4+   2,ra 


hr)(n+    1)-A         (n-    1,4+   2,m  +    2 


V-  (4+    1)  -A  (n+    1,4  ,m+    2 

A  (r>+    1  ,4+   2  ,m+    1 

s  ,  c 


^ra  •  B         (n+l,4+3,m 


s  ,c 


?/*ra-B         (n+    2,4+    2,ra 


mn  •  B         (  n  ,  4  +    3  » m 
s  ,c 


m4  •  B         (n+    2  ,4+    1  ,m 
s  ,c  v 


(A-6) 


where    it    is    assumed  A    ,     B     are    chosen    when    the    left    hand   side    is 

s         s 

<S|  (H-E)  l\  >  and   A    ,     B      when    it    is    <c|  ( H-E)  |  \>. 

3 •       Determination    of    the    Required    Integrals 

Clearly    not    every    term    in    the   above   matrix    elements    requires    a 
different    integral.      Many    of    them   are    identical    or    have    zero 
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coefficient.   For  example,  for  the  34  element  basis  set  (M  =  5), 
while  there  are  595  independent  elements  each  in  H  and  S,  repre- 
senting 52  terms  and  4  terms  each,  respectively,  or  a  total  of 
more  than  33,000  terms,  only  about  700,  or  slightly  more  than  two 
percent,  use  different  integrals  and  all  700  different  integrals 
can  be  generated  using  the  recursion  relations  proved  in  Appendix  B 
by  actually  carrying  out  the  evaluation  of  about  100  integrals  for 
which  /i  =  2.   Similarly,  in  evaluat ing  <S |  ( H-E) | x  > and  <c] H-E) | x  >  . 
where  there  are  34  elements  in  each  expression,  each  element  in  turn 
made  up  of  26  terms,  for  a  total  of  1768  integrals  required,  only 
268,  or  slightly  less  than  one  sixth  are  distinct  and  again,  using 
the  recursion  relations  only  103  terms  for  jj,    =  2  need  be  evaluated 
explicitly  and  they  can  all  be  found  from  3°  different  numerical 
integrals.   Nevertheless  it  is  the  numerical  evaluation  of  these 
integrals  which  overwhelmingly  dominates  the  time  required  to  solve 
the  complete  problem  and  hence  the  requirement  only  to  evaluate 
those  terms  needed. 

To  insure  that  no  terms  were  missed  in  computing  the  integral 
tables,  preliminary  computer  programs  were  written  to  examine  each 
of  the  matrix  elements  and  to  display  in  tabular  form  which  integrals 
were  needed  to  construct  the  matrix  elements  for  a  given  basis  set 
size.   Tabulated  were  all  distinct  integrals  required  for  which  the 
coefficient  was  non-zero  in  at  least  one  case  plus  those  required 
in  order  that  the  recursion  relations  could  be  used  to  generate 
the  terms  required  for  \i    =   2 .   From  this  information  the  integral 
table  programs  were  constructed  so  that  only  those  terms  required 
were  actually  calculated,  and  none  that  would  be  required  was  omitted. 
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APPENDIX  B 
Integral  Recursion  Relations 
1.   Recursion  Relations  for  u  >  2[5l] 


and 


Consider  the  integrals 

A(V,X,H)=  5d"ld^2ri"2r2"2ri22F(ri'r2)  (B_1) 


B(V,A,H)  =   \d^1dVl"2r2"2ri22F(riA:2)COS912  (B'2) 


and   rewrite    (B-l)    in    the    form 

using   the    law   of   cosines.      This    transforms    immediately    to 

A(V,X,U)  =    A(V+2,X»M--2)+    A(V,X+2,u-2)  -    2B(  V  +  l  ,X  +1  ,H  -2  )  .  (B-3) 

Similarly    B(v,X>^)    becomes 

B(v,X,u)  =    B(V+2,X,H-2)  +    B(v,X+2,u-2)-   2\  r^_1r  *"  V^~  V(  r  x  ,  r  2)  cosQ^dr  .^h:  2< 

Let  the  integral  in  the  third  term  be  designated  I,  then 

I  =  A(Vl,X*l^-2)-5r^-1^-:t1r^-V(r1,r2)sin2ei2<te1d?2. 

Now  recall  that  the  angular  part  of  the  dr  dr   integration  is 

simply  sinO   d9   ,  whence  the  angular  part  if  the  second  term  in  I  is 

S(V  r2"  2lir2COS912)^"  Sin2912d912 

=  ^2  -TT  f  S-  «1+  V  2rlr2COS912)^  ]si"2ei2d0i2 
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which  may  be  integrated  once  by  parts  to  give 


Tr-S 


r^~   cos©, „sin9,  d9. 


U--2  r  r   J  12       12     12   12 


and   therefore 


I    =    A(V+1,X+1,ll-2)+    JL    B(V,A,^) 

Li  -z 

from  which  it  follows  immediately  that 

B(v,X,'j)  =  —^   |b(v+2,X  ,u-2)+B(v,A+2,M.-2)-2A(v  +  1,X+1^-2)}       (  B-4) 

2 .       Recursion    Relations    for    ii    =    1 . 
When  u    =    1,     (B-l)    becomes 

C    \)-2  X  -2  1         -      - 

A(vA,l)=  yl      r2      F(r1,r2)    j^-   dr^    . 

Now  recalling  that  in  the  coordinate  system  in  use  r  >    r  ,  so 
1/r    can  be  expanded  in  a  series  of  Legendre  Polynomials 

I 

co       y 

-^—  =  £    -7TT  Pfl(cosO,_) 

ri2    x=°  r"+1   £      12 


and  the  angular  integration  becomes 

I  t           ,1 

-    r     fn  -    r^     1 

£    -f-r  \  P,(cos9,-)sin9._d9._  =  S    -*-=-  \   P,(u)du. 


But    P    (u)    =    1    so   the    integral    on    the   right    can    be    rewritten 


S   ,    P0(«)Pi(»)du 


-1 

and  by  virture  of  the  orthogonality  of  the  Legendre  polynomials, 
all  terms  in  this  expansion  vanish  except  for  I    -   0 .   Therefore 

A(V,X  ,1)  =  A(V-1,V,2) .  (B-5) 
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In  a  similar  expansion  in  the  case  of  B(v,X,l)  all  terms  but 
I    =    1  vanish,  giving 


B(V,X  ,1)  =  |  A(V-2,X+1,2)  . 


(B-6) 


3.   Recursion  Relations  for  A    (v,X,ll)  and  B    (v,X,u) 

Referring  to  (53a) 

As(V,X,^)  =  A1(v,X»M-)+  A2(V,X^) 

and   applying    (B-3)    above 

A    (V,X»U)=   A1(X+2,v,^-2)+   A  (X,V-2,p.-2)-    2B    (X  +1  ,  V  +  l  ,M-  -2  ) 

+   A2(V,X+2,u-2)+    A2(V+2,X  ,M--2)-    2B2(V+1  ,X+l,M-2) 

=    A    (\)  ,\+2,u-2)  +   A  (V+2,X,U-2)-    2B    ( V+l  ,X +1  ,M<-2) 


and   similarly    with    B   (v,X>^)-      Thus    A  and    B  satisfy    the 

s v  '  s  ,c  s  ,c 

normal    recursion    relationships    for    U-    >   2  .       For    \i    =    1 

As(v,X,l)    =  A1(X,vM)+  A2(v,X,l) 

=    A1(X-1,V,2)+   A2(V-1,X  ,2) 


but 


so 


A2(V-1,X,2)    =    As(v-l,X,2)-    A2(X,V-1,2) 


A         (V,X,1)=AC        (V-1,X,2)+A  (X-1,V,2)-    A  (X,V-1,2)  (B-7) 

where   A,     is    chosen   when   A      is    to   be    evaluated   and   A„    is    chosen   when 
Is  3 

A      is    to   be    evaluated.       For    B        (v,X>±)    the   same  procedure   yields 

Bs,c(V'X,1)=    3   {As,c(;~2'X+1'2)+    A1?3(*-2>V+1,2)-   A1    3(X+l,V-2,2)} 

(B-8) 
which  completes  the  derivation  of  all  the  relations  needed  to 
evaluate  the  matrix  elements  listed  in  Appendix  A. 
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APPENDIX  C 

Rayleigh-Ritz  Variational  Method 

Let  an  arbitrary  quadratically  integrable  wave  function  Y  be 
expanded  in  terms  of  the  eigenfunct ions  of  a  bound  system 

Y  =  2  aX-  (C-l) 

1  x    x 

where 

(H-E)|\.>  =  0 

then  on  the  function  ¥  the  expectation  value  of  H  is 

<Y|h|Y>  =  E    a*a  .<X-|h|\  > 
i,j  x   3      x     J 

=  ^U.|2E.     .  (C-2) 

If    E      is    the   oround   state    of   the    system,    then    E      ^    E.    so 

o  O  1 

E  Ela.l2    *   S  !a.|2E. 


Ol  1  ill 


whence 


assuming 


But 


<  y|h|y  >  *   E 


<  YlY  >  =  1 


<  y|h|y  >  =  E^ 

the   expectation    of    H   on   Y  ,    so    if   Y    is    found   so  as    to   minimize   E~, 
then   Ey   ^    E      and   E,^    provides    an    upper    bound   on    the   value    of   E    . 
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In   general 

By   =  <   y|h|y   >/<Y|Y>  .  (C-3) 

If   now   a    trial    wave   function   Y       is    chosen   and    expanded    on   a    suitable 


basis   r\  . 
1 


whence 


N 


Y      =  E        b.rj  .   ►  Y  (C-4) 

S        b*b  <T7 .  |  H 1 77     >  £        b*b.H.. 

.       .         l     J        *  J  .       .         1     J     11 

p       =    izJ =  U 


t  * 

E         b.b  .S.  . 
i,j       *    J    iJ 


E        b .  b  .<ti  .   n  . 
i,j      i  J    'i'  fJ 


where   H.  .    =  <  77  .  I  h|t?  .>     and   S..    =  <T7  .  |"»7  .>  .       Now    choosing    the  b's 
to    satisfy 

5bT  =  °>         i  =  L....N 

l 
gives 

■5—   E        b.b  .S.  .+    E      s-r—  E        b.b.S.  .    =   ?—-  E         b .  b  .H .  . 
bb.    io-      x   j   ij        t    bb.    ifj      i    j    ij        6b.    ifJ      i   j    ij 

from   which    arise    the    N    linear    equations 

E(H.  .-    E.S.  .)b.     =0,  j    =    1,...,N  (C-5) 

i       ij         t    ij '     l 

which  is  just  the  first  part  of  the  Harris  expansion  technique, 

wherein  the  smallest  eioenvalue  E   gives  an  upper  bound  on  E  . 

t  '  '  o 

The  above  method  extends  readily  to  the  excited  states  [ 52 J . 
This  is  achieved  by  choosing  a  trial  wave  function  which  is 
orthogonal  to  all  the  states  lying  below  the  desired  one.   Such  a 

function  must  take  the  form 

N-l 
Y  =  cp  -  E    X  <X  I cp  >  (C-6) 


n 


=0   n   n 
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where  cp    is    arbitrary    and   the  \      are    the   lower    lying    eigenstates. 

Then    expanding  Y    according   to    (C-l)    it    is    seen    that    the   a.    are 

zero   for    i    less    than    N,    the    index   of   the   state   desired,    from   which 

it    follows    that    E,.   ^    E„  . 

t     N 

Hylleraas  and  Undhcim  [32]  have  shown  that  this  condition  is 
automatically  met  by  (C-5),  the  n    eigenvalue  of  that  equation 
being  greater  than  or  equal  to  the  energy  of  the  n    level  of  the 
system.   Thus  it  is  seen  that  an  important  bonus  in  the  Harris 
expansion,  method  is  a  cataloging  of  the  energy  levels  of  the  bound 
system  formed  by  the  scattering  atom  and  the  incident  electron. 
Included  in  this  cataloging  are  the  autoionizing  levels  of  the 
bound  system  which  give  rise  to  many  of  the  observed  scattering 
resonances  [44J  • 
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APPENDIX  D 


Discussion  of  Numerical  Calculations 


1.   The  Eioenvalue  Problem 

From  Section  II  C.l,  a  solution  for  eigenvectors  and  eigen- 
values is  sought  for  the  equation 

(H  -  \S)c  =0  (D-l) 

Since  S  is  not  the  identity  matrix  (D-l)  is  not  in  the  form  usually 
associated  with  eigenvalue  problems  occurring  in  physics  where  the 
space  is  orthonormal  and  S  reduces  to  the  identity  matrix.   This 
generalized  eigenvalue  problem  has  been  extensively  discussed  in 
the  mathematical  literature  L  4l D -   In  the  present  case,  however, 
it  is  possible  to  reduce  the  problem  to  standard  form  and  treat 
it  by  the  simpler  methods  available. 

Consider  first  the  matrix  S_.   Its  elements  are  the  inner  pro- 
ducts of  the  chosen  basis  vectors 

sij  =  <xiIxJ->    •  (D~2) 

Since  the  basis  functions  are  real  S.  .  =  S..  and  therefore  there 
exists  a  unitary  transformation  which  diagonalizes  S  L53]>   Physi- 
cally this  is  equivalent  to  transforming  to  a  new,  orthogonal  set 
of  basis  vectors.   In  this  new  set 

S!  .  =  5.  .  <X!|X'> 
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3 
and   S ! •>    0    for    all    i.  But    diagonal iz ing   a    matrix    in    this    manner 

11 

means    that    the    new  matrix   has    as    its    diagonal    elements    the    eigen- 
values   of   S.      Therefore    it    can   be   concluded    that    all    the   eigen- 
values   of   S   are  positive   and   hence   S_  is   positive   definite. 

Now  consider  the  positive  definite  quadratic  form  associated 
with  S,  x  -S/x.  There  exists  a  real  non-singular  transformation 
U  and  a    vector    y    such    that   x    =   U       *y    and   such    that    (Ref .    53,    p.    337) 

— »T         -T  -1      -*  ~T     — » 

y    -U      -S-U      .y    =   y    -y  (D-3) 

"T    "•  -T      -1  T 

called  the  "cannonical  form"  of  x  -S'x,  where  U    =  (U   )  .   Hence 


u^.s-u'1  =  x  (°-4) 

where  I  is  the  identity  matrix.   So 

S  =  uTu  .  (D-5) 

In    other    words,    a   positive    definite    symmetric    matrix   can   always    be 

factored    into    the  product    of   a  real    non-singular    matrix  and    its 

transpose.       But    because   S    is    symmetric,    this    factorisation    is    not 

unique,    for 

N 

s. .  =  s      u.  .u.  . 

lj         k=1      ki    kj 

2 
represents    n      equations,    only   n(n+l)/2    of  which  are    independent. 

Therefore   one    is    free   to   choose   n(n-l)/2    of   the   U^      arbitrarily, 

*j  m 


3 

Unless    one    of    the   new   basis    vectors    is    null,    which    would   mean 

the    dimensionality    of    the    space    is    reduced    -   a    possibility    which 

can    be   excluded    from   the  present   problem. 
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and    it    is    convenient    to    choose 

U,       =0        JL    >    m 

Jem 

for    then   U   becomes    an   upper    triangular    nuitrix    (one    in   which   all 
elements    below    the  principal    diagonal    are   zero). 

With    the   above    choice   of   factorization    of  S_,     (D-l)    can   be 
rewr  it  ten 

U_T-(H-  ^jJT*U)U~    -U-c    =   0 
and    the    original    problem   reduces    to 

(if^H-lf1  -  \l)-x    =0  (D-6) 

~*  ~*  .  .....  -T     -1 

where  x  =  U  •  c  .   Since  H  is  symmetric  it  is  trivial  to  show  U   -H-U 

is  also  symmetric,  for 

-T     -1  -T     -1         -1     -T 
(U   -H-U   )   .  =S    U..H..U..   =E    U.  H.  .U., 
'mj    i)k   ik  ij    ik    ijk   im  ki  jk 

=  £    U~*H.  .Ul1  =  (U"T-1!'U_1)  . 
i,k   jk  kj       v- jm 

and  the  problem  has  been  reduced  to  the  standard  problem  of  a  real 
symmetric  matrix  in  an  orthonorraal  space. 

The  advantage  of  the  form  of  the  factorization  of  S  chosen  is 
that  because  U  is  triangular  its  inverse  is  particularly  easy  to 
obtain.   The  factorization  is  done  by  the  square  root  method  of 
Cholesky  L54l  and  U    is  then  found  by  solving  the  algebraic 

equations  [55] 

U?*U..    =    1  1 

k         ""  V         1-1.....N 

e   u:V  =o,  i<kj  (D-7> 

for    the   u7 .    . 
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Because  the  range  of  values  of  the  elements  of  S^  and  H  was 
often  several  orders  of  magnitude  (sometimes  as  many  as  fifteen 
or  more)  an  additional  check  on  the  accuracy  of  U    as  found  by 
the  above  method  was  used.   If  the  U    found  above  is  taken  as  a 
first  approximation  to  the  exact  inverse,  then  the  following 
iteration  scheme  can  be  used  to  improve  the  value  of  the  inverse  L  56  ]  . 

Let  C.  be  the  i    approximation  to  the  exact  inverse,  U   ,  and 
define  the  error  matrix 


Then 


R.  =  I  -  C. -U  .  (D-8) 


C.  .  =  C.+  R.C.  =  C.+  (I-  C.   •  U)-C. 
—l+l    —  l   —l—i    —l    —  —l    — '  — 1 


and 


R.  .  =  I  -  [C.+  (I  -  C.-U)-C.}'U 

—1+1    —     —l    —  —  l  — '  —  l   — 

=  (I"  £i-I')2  =  £                                                 <D"Q) 

v  —  r>' 


wher  e 


R   =  I  -  C  -U 
— o    —   — o  — 

and  C   is  taken  to  be  the  value  of  U    found  by  the  Cholesky 
— o  — 

method  above.   Clearly,  if  the  norm  of  R   is  less  than  one  the 

— o 

method  converges  quite  rapidly,  and  since  the  value  of  C   is 

already  quite  accurate,  ||R  ||«  1.   In  practice  no  more  than  one 

-20 
iteration  of  C   was  ever  required  to  make  the  elements  of  R .<  10 
— o  M  1 

The  eigenvalue  problem  was  solved  using  the  Jacobi  variable 
threshold  method  L4l>55j-   This  method  proved  to  give  the  most 
accurate  results  although  it  is  known  to  be  considerably  slower 
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than  other  methods  tried.   The  accuracy  of  the  eigenvalues  was 
verified  by  the  following  method.   For  a  given  eigenvalue,  X  ,  the 

r- 

determinant  |  H-X  S]  was  formed  and  compared  with  the  determinant 
r 

when  X   was  replaced  by  X  +  5  where  6  was  typically  10    (compared 

to  X   which  was  on  the  order  of  one).   Then  X   was  considered  cor- 
[1  U 

rect  to  six  places  if  the  respective  determinants  were  related  by 

|H-(Xn  +  6)S|<|H  -  X  s|<|H  -  V.  +  6)s|  (D-10) 

and  the  first  and  last  terms  above  were  of  opposite  sign.   Accuracy 
of  the  eigenvectors  was  not  so  clearly  established,  though  a 
similar  procedure  to  the  above  was  followed.   The  product  (H-Xs)c 
was  formed  and  it  was  noted  that  as  X  was  perturbed  as  above,  all 
elements  in  the  resultant  vector  changed  sign.   While  this  indicated 
that  the  eigenvectors  probably  are  reasonably  accurate  no  quanti- 
tative bound  could  be  placed  on  their  error.   Further  indication  of 
the  accuracy  of  the  eigenvectors  can  be  inferred  from  the  close 
agreement  of  the  phase  shift  results  for  hydrogen  with  the  accepted 
values  found  by  Schwartz  [l2]. 
2 .   Numerical  Integration 

From  the  beginning  it  was  felt  that  the  simplest  and  most 
reliably  accurate  means  of  evaluating  the  necessary  integrals 
numerically  would  be  some  form  of  Gaussian  quadrature  (7 5 7 J  -   Any 
numerical  integration  scheme  attempts  to  replace  the  integration 
by  a  finite  sum  so  chosen  as  to  make  the  error  tolerably  small 


s 


n 
f (x)dx  «  S    w.f(x.)  (D-ll) 

i=0   1    x 
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where  the  w.  are  weights  to  be  assigned  to  the  function  evaluated 
1 

at  the  sampling  points  x.,  in  other  words  a  weighted  average  of  the 

function  at  the  various  sample  points.   If  the  x.  are  fixed,  equally 

spaced  points  then  the  n  +  l   w .  can  be  chosen  to  define  a  polynomial 

of  order  n  with  which  to  approximate  f (x) ,  and  if  f (x)  is  itself  a 

polynomial  of  order  n  or  less,  the  integration  will  be  exact.   If 

however,  one  is  free  to  vary  the  x.  also,  there  will  now  be  2n  +  2 

1 

parameters    available   and    it    will    be  possible    to   define   a   polynomial 

of   degree   2n    +   1    with   which    to   approximate    f(x).      Thus,    given    the 

values    of    x.    and   w.    in    each   case,    the   same  amount    of   labor    results 
i  l 

in   a   more   accurate    integration    in    the   second    case. 

Adopting    this   approach,    one    can   approximate   f(x)    with   a 
Lagrangian    interpolating  polynomial    of   degree    (n    +    1)    of    the   form 


f(x)     =    g(x)F(x)=   S         L    (x)g(x)F(x    )+    R     (x)  (D"12) 

1=0       r  in 

where 

n 

L.(x)    =7T     (X"XJ}  (D-13) 

1  j=0    (x.-x.) 

and  R  (x)  is  the  remainder  term,  given  by 

n 

R(x)  =7T  (x-x  )F(n  +  1)(Ug(x)/(n+l)  !,      a  <§>  b        (D-l4) 
i=0 

The  reason  for  factoring  f(x)  into  g(x)F(x)  will  become  apparent 

below.   Now  if  F(x)  is  a  polynomial  of  degree  2n  +  1  then  F        (§) 

must  be  a  polynomial  of  degree  n.   Letting 

F(n+1)ri 

-7 — rrt^  =  q   (*) 

(n+1 )  I  irtK     ' 
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where  q  (x)  is  a  polynomial  of  degree  n ,  (D-l4)  becomes 


n 


Rn(x)  =  qn(x)ff      (x-x.)g(x)  (D-15) 

i=0 


Integrating  (D-12)  gives 


pb  n        >-«b  pb  n 

\  g(x)F(x)dx=  £   F(x.)\  g(x)L  (x)dx  +  \  g(x)q    (x)ff   (x-x.)dx    (D-l6) 

Ja  1=0    X  Ja  Ja      n    i=0 


Now  (D-l6)  is  of  the  form  (D-ll)  where 

#-«b      n   ( x-x  . ) 

vvi  =\  9(x)7T  ^-r^  dx  (D-i7) 

a   j=0   ^  r 
jVi 

except  that  the  x.  have  not  been  chosen.   The  object  is  to  so 

choose  the  x.  that  the  second  term  of  (D-l6)  vanishes.   This  is 

n 
easily  done  if  q  (x)  and  J~f     (x-x.)  are  eqoanded  in  a  series  of 

n         i  =0       J 
polynomials  orthogonal  on  [a,b]  with  rcspoct  to  the  weight 

function   g(x) .   A  weight  function  appropriate  to  a  particular 

set  of  orthogonal  polynomials  is  often  a  natural  factor  in  the 

function  to  be  integrated  and  therefore  the  particular  expansion 

to  be  chosen  may  be  dictated  by  the  problem,  hence  the  utility  in 

carrying  the  factor  g(x)  through  the  derivation. 

Let  the  polynomials  chosen  for  the  expansion  be  designated 

P  (x)  whereupon 

7T  <*-* i>  =H  Wx)  <D"l8) 

i=0         x_u 


and 


q  (x)  =  E    c  P  (x)  (D-19) 

n       m=0   m  mv  ' 
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Combining  these  two  expansions  in  (D-15)  and  applying  the  ortho- 
gonality condition 

.b 


{    g(x)Pjt(x)Pm(x)dx  =  0,  I    jt   m  (D-20) 


reduces    the    error    term    to 


\    g(x)qn(x)ff      (X-X    )dx    =S      be      \     g(x)[P    (x)j     dx 

Now  if  the  x.  are  chosen  to  be  the  zeroes  of  P   _ (x) ,  then  the 

1  n  +  1       ' 

expansion    (D-l8)    reduces    to 

77"    (x-x.)    =   b        P      .  (x) 
l\      v  \>  n+l    n+lv     ' 

i=0 

i.e.,  b.  =  0   O^j^n,  whence  (D-21)  vanishes,  and  the  integration 

is  exact.   Thus  evaluating  the  function  at  n+l  carefully  chosen 

points  suffices   to  integrate  it  exactly  if  it  can  be  expressed 

in  terms  of  a  polynomial  of  degree  2n  +  1  or  less,  a  considerable 

improvement  over  the  equally  spaced  interval  methods,  as  was 

promised  above.   The  one  remaining  problem  is,  then,  to  evaluate 

the  w.'s  and  find  the  x.'s.   Fortunately,  this  has  been  done  for  a 
x  1 

large  class  of  orthogonal  polynomials  by  Stroud  and  Secrest  [58] . 

Although  the  integration  limits  and  the  form  of  the  integrands 
( 49) >  (  5l)  >  (56)  and  (57)  suggest  the  use  of  the  Laguerre  poly- 
nomials for  the  numerical  integrations,  the  almost  periodic  nature 

of  the  functions  and  the  ever  increasing  amplitude  in  the  absence 

—  x 

of  the  e  '  factor  dictated  the  use  of  a  finite  interval  formula 

where  several  intervals  could  be  integrated  separately  and  then 

—x 

added  together  to  give  the  final  result.   Since  the  factor  e 

forces  the  integral  to  zero  at  large  x,  examination  of  the  integrands 
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indicated  that  an  upper  limit  on  the  integration  of  x  =  100  would 
introduce  a  negligible  error  from  the  neglected  region  in  the  worst 
possible  case.   Because  the  variable  change  to  conform  to  the  con- 
vergence limits  of  the  polynomials  used  in  the  expansions  (D-l8,19) 
is  simple  and  because  the  weight  function  is  unity  the  decision 
was  made  to  use  Legendre  polynomials  in  the  integration. 

Because  of  the  almost  periodic  nature  of  the  integrands,  when 
the  instantaneous  frequency  of  1he  sinusoidal  variation  is  large  the 
areas  under  the  integrand  above  and  below  the  axis  are  nearly  equal 
and  substantial  loss  of  significance  occurs  when  these  two  values 
are  subtracted.   Often  as  many  as  eight  to  twelve  significant 
figures  can  be  lost  due  to  subtraction,  hence  great  care  must  be 
taken  to  avoid  the  problem.   The  approach  adopted  was  to  seek  a 
function  whose  exact  integral  is  known  and  which  approximates  the 
desired  integral  closely  enough  so  that  the  difference  between  them 
is  small,  then  evaluating  this  small  difference  numerically  and 
adding  the  result  to  the  known  value  of  the  approximating  function. 
This  method  can  be  compared  to  the  measurement  of  the  difference 
between  two  large  quantities,  say  frequencies.   Often  the  difference 
can  be  determined  quite  accurately  even  though  the  absolute  magni- 
tude is  known  only  approximately. 

In  executing  this  approach  to  the  integration  of  the  functions 
at  hand  the  first  step  is  to  find  all  the  zeroes  of  the  integrand 
between  arbitrarily  chosen  minimum  and  maximum  values  of  x.   These 
values  are  then  used  as  the  end  points  of  the  subintervals  in  an 


integration  of  the  form  (D-22) 

Zl         n     Zi+1  Z-+i  10° 

I(x)=^   f(x)dx+  £  =  |\     [f(x)-hi(x)Jdx+\     hi(x)dxj+V    f(x)dx 


z  .  ~z  .  z 

i  i  n  +  1 
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where  the  z.  are  the  n  +  1  zeroes  of  f(x)  found  above  and  the  h.fx) 
1  '  x      ' 

are  functions  whose  exact  integral  is  known  and  which  approximate 
f(x)  over  the  interval  indicated.   The  h.(x)  are  chosen  as  indicated 
in  Figure  23.    Note  that  h.(x)  is  amplitude  modulated  from  segment 
to  segment  so  that  f(x)  -  h.(x)  is  always  mostly  positive  and 
smaller  than  f(x)  by  at  least  an  order  of  magnitude  regardless  of 
the  sign  of  f(x).   It  is  this  difference  function  which  is  inte- 
grated numerically,  and  since  h.(x)  is  a  different  function  for 
each  half  cycle  of  f(x)  the  numerical  integration  is  done  in 
segments  between  zeroes  of  the  function  using  a  32-point  Gauss- 
Legendre  quadrature  on  each  segment. 

The  form  chosen  for  h. (x)  was 

l 

h.(x)  =  (1  +^-  )sin(Ax+B)e'Xxy(l-e"PX)3  (D-23) 


where  p  is  +1  if  f(x)  is  negative  and  -1  if  f(x)  is  positive, 

B  =  1  if  I   or  I   is  to  be  found  and  B  =  0!/2  ( Z+a)  if  I,  or  T„  is 
r         s     c  r  1     3 

to  be  found,  and 

i+I      i 
so   chosen      that         sin(Ax    +    B)    vanishes    at    the    end  points    of   each 
segment   and   nowhere    in   between. 

When    chosen    in    this    manner    the    integral    of   h.(x)    alone    is 

-Bx    3 
exact   and   has    a    rather    simple   form.       If    (1    -    e    *"     )       is    expanded 

four    integrals    of   the   form 

P    i+I  ~a/jx 

S£    =   \  sin(Ax+    B)e  x7    dx ,    i    =    1,2,3,4)  (D-25) 

z  . 

l 
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f(x)-h.+1(x) 


i+2 


/hi+i<*> 


Figure    23. 
Example    of   Method   of   Choosing   h.(x) 


result,    where   a.     is    1  ,    1    +    p  ,    1    +   2p    and    1    +    3p    in    turn 
integrate    immediately    to 

y+l 


These 


Si    = 


j=l    (-D)J(y+i-j)!    V  1+1  y  (D'26) 


where   9   and  p    satisfy    the   relationship 

-a    +    (~1)2A   =  p(cos0    +    (-l)2sin9) 
and    finally 


Zi  +  1 


^  h.(x)dx     ::     P(1+X-)(S1-3S2+     3S3-     S4) 


(D-27) 
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Express  ions  (D-26 ,27)    are   valid    for    y^  0.       During   the    course    of 
evaluating    the  phase   shifts,   however , it    is    necessary    to   evaluate 
the    various    integrals    for    y   -   -1 .       For    this    case   a    slightly    dif- 
ferent  scheme   must    be   used    to   evaluate   h.(x).       If   one   replaces 

-Bx 
(1    -    e   p'  )/x    in    (D-23)    by    the    integral    identity 

-Dx  rl 

1    -    e  h       =    B    {    e"'Xydy  (D-28) 

then 

h.(x,    =  p(1    +  I '  ,   C'sintAx.Bje^d-e-P^^-P^dy  <D"29> 

i  iU    Jo 

Carrying   out    the   x   integration   first    and  making    the   appropriate 
variable    changes    yields 


B(J&+1)+1      -z.w        -y- ■  ^ 

rv  '  l  r+1 


w 


S 


■S 


,    i-^ dw  (D-30) 


6^+1  w      +    A 

and    finally 
z  •     -i 

r  1+1  p 

\  h.(x)dx    =    -    P(l     +    — )-A'(S    -2S    +S0) 

J  iV     ;  V  10;  o         1       2'  (D-31) 

'i 
Unfortunately    (D-30)    must    be    integrated   numerically    also,    however 
it    is    a    relatively    slowly    varying    function    of   w  and    is    always 
positive    so    integration    by   a    32-point   Gauss -Legendre   quadrature 
gives    more    than   adequate   accuracy. 

To   complete    the    integration    of    (D-22)    the   first   and    last    terms 
were   evaluated   using         single    512-point   Gauss -Legendre    quadratures. 

Since    the   accuracy    of   the   numerical    integration    is    crucial    to 
the   accuracy    of    the  phase    shifts    calculated   for    helium    ion    scat- 
tering,   a    great    deal    of    effort   was    expended    to    insure    that    the 
methods    were   both    correct   and   internally   accurate. 
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The  first  and  most  obvious  check  is  to  increase  the  number  of 
sample  points  in  the  numerical  integrations.   This  was  done  by 
increasing  the  quadrature  method  on  each  interval  to  64  points. 
Changes  to  the  integrals  on  doing  this  were  typically  in  the  12 
to  14    digit,  indicating  that  the  differences  were  principally 
due  to  round-off  error  since  the  maximum  number  of  digits  available 
was  16 . 

The  integration  scheme  described  above  grew  out  of  a  similar 
method  which  was  used  extensively  prior  to  the  present  one  for 
the  same  calculations.   The  present  method  was  adopted  in  an  effort 
to  extend  the  useable  range  of  the  parameters  involved  in  the 
integrands.   The  problems  with  the  former  method  arose  principally 
because  it  was  not  sufficiently  accurate  when  the  oscillation  fre- 
quency of  the  sinusoidal  function  was  high.   However  in  the  region 
where  both  schemes  were  adequate  (which  proved  to  be  almost  every- 
where) the  integrals  were  consistently  r epr oduceable  to  12  to  14 
places.   The  choice  of  integration  intervals  and  of  the  h.(x)  in 
the  earlier  method  was  sufficently  different  that  it  constitutes 
an  almost  independent  check  on  the  internal  accuracy  of  the 
integrals . 

Three  different  methods  were  used  to  make  the  integrals  exactly 
integrable  so  that  the  numerical  results  could  be  compared  directly 
with  known  answers  and  to  provide  a  check  on  the  correctness  of  the 
expressions.   The  first  was  to  let  k  -*  °°  while  holding  the  ratio 
k/a  constant.   This  forces  the  logarithmic  terms  to  zero  and  re- 
duces the  integral  to  a  standard  form  which  can  be  integrated 
directly.   The  direct  evaluation  was  done  on  a  Hewlett-Packard 
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Model  9S10A  computer  which  gave  results  to  ten  significant  figures. 
These  results  were  then  compared  with  the  numerical  results  and 
except  in  the  cases  where  k/a  was  at  its  extreme  limits,  agreement 
was  obtained  to  eight  to  ten  places.   This  method  was  tried  only 
with  the  older  integration  scheme  and  so  some  decay  at  the  limits 
of  k/0!  was  expected.   In  addition,  the  transition  from  the  form 
of  the  integral  when  k  was  finite  to  its  form  as  k  ~*  °°  appeared 
smooth  and  with  no  obvious  discontinuities  or  other  strange  behavior 
in  the  transition  region. 

A  second  check  took  advantage  of  the  fact  the  computer  program 
was  written  to  allow  for  variation  of  the  nuclear  charge,  Z. 
Although  the  intent  was  to  limit  the  value  of  Z  to  integers, 
nothing  in  the  program  or  the  mathematics  forbade  varying  the 
charge  continuously  from  one  to  two.   When  this  was  done,  it 
was  again  seen  that  the  result  for  Z  =  2  transitioned  smoothly 
into  that  for  Z  =  1  with  no  discontinuities  as  the  charge  was 
var  ied . 

Since  for  Z  =  1  the  integrals  all  become  exactly  integrable, 
the  computer  program  did  not  do  the  calculations  required  for  the 
numerical  integration  in  the  case  of  Z  =  1 ,  but  skipped  to  a  sub- 
routine which  evaluated  directly  the  exact  integrals.   The  final 
comparison  method  was  to  circumvent  the  exact  integration  routine 
for  z  =  1  and  do  these  integrals  numerically,  using  the  method 
described  at  length  above.   Again  the  agreement  between  the  two 
calculations  was  excellent,  the  integrals  and  the  phase  shifts 
typically  being  identical  in  the  first  twelve  to  fifteen  places. 
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The  final  check  was  to  test  the  sensitivity  of  the  entire 
calculation  to  perturbations  in  k.   It  was  noted  that  even  at  the 
extreme  limits  of  k/Q!  where  accuracy  was  expected  to  be  poorest, 
only  some  of  the  integrals  showed  evidence  of  toss  of  significance 
from  the  subtraction  problem  mentioned  above.   For  these  integrals, 
however,  a  small  perturbation  in  the  value  of  k  may  have  a  rela- 
tively large  effect  on  the  integration  and  if  the  inaccurate 
integrals  dominated  the  computation  then  they  could  significantly 
distort  the  results.   As  before  when  this  perturbation  was  arti- 
ficially fed  into  the  problem  its  effect  was  well  within  the  limits 
considered  acceptable  (results  unchanged  to  six  places  or  more  for 
a  perturbation  of  k  in  the  seventh  place)  except  at  the  extremes. 

Since  none  of  these  tests  indicate  any  serious  problems  with 
accuracy  it  is  felt  that  the  errors  inherent  in  numerical  inte- 
gration have  been  kept  well  within  tolerable  bounds  in  this  problem 
and  the  results  can  be  relied  upon  with  considerable  confidence. 
It  also  proved  to  be  true  that  those  regions  where  accuracy  was 
lost  were  easily  avoidable. 
3  •   Evaluation  of  V ( 1  -  ia) 

From  the  definition  of  the  gamma  function 

CO 

T(P)  =  \  e'XxP_1dx  (D-32) 

J0 

one  finds 

00 

T(l-ia)  =  \  e"Xx"iadx  (D-33) 

J0 

but  e   *  '=  x,  so  (D-33)  can  be  rewritten 

r(l-ia)  =  \  e"Xcos(a^Jx)dx-i  \  e~Xsin  (a^  x)  dx  (D-34) 

=  A  +  iB 
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whence 


-1  B 


rg  [r(i-ia)}  =  tan'   -  (D-35) 


A 


and  the  evaluation  of  the  complex  gamma  function  is  reduced  to 
evaluating  two  real  integrals.   The  integrals  have  to  be  evaluated 
numerically,  however,  and  are  not  suitable  for  this  in  the  form 
(D-34)  because  the  rapid  oscillations  of  the  integrand  near  the 
origin  could  contribute  to  a  substantial  error  in  the  numerical 
integration.   To  overcome  this  problem  note  that  an  integration 
by  parts  of  each  of  the  integrals  (D-34)  introduces  a  factor  x 
into  the  integrand  and  this  has  the  desirable  effect  of  reducing 
the  contribution  of  the  integrand  near  the  origin.   Each  re- 
petition of  the  partial  integration  introduces  an  additional 
factor  of  x  into  the  integrand.   Carrying  out  this  process  three 
times  transforms  the  integrals  (D-34),  after  considerable  algebra, 
i  nto 


B  =  - 


(l+a2)(4+a2)(9+a")^0 


and 
A  = 


(l+a2)(4+a2)(9+a2),J0 


r°°  3  -x  /6   2        ,      2  1 

-  \  x' e   1  —  (a  -l)sin(avjc)-(a  -11 )  cos  (a.":x)  f  dx 

(D-36) 

V  x3e~X|(a2-ll)sin(a>.x)+  ~(a2-l  )cos  ( a."".x)}dx 

(D-37) 


which    forms    are    now   quite   suitable    for    numerical    integration. 

The    integrations    of    these   functions    were   carried    out    using   a 
256-point   Gauss -Legendre    quadrature    over    15  arbitrarily    chosen 
intervals.      To    check  accuracy    the   results    were    compared   with 
another    program    which    calculated   the    complex   gamma    function    using 
a    Pade-power    approximation    method  and    claimed   9-place   accuracy. 
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The    internal    accuracy    of    the    integration    scheme   was    checked   by 
comparing    the   256-point    integrations    to    512-point    integrations 
on    the    same    intervals.      This    comparison    indicated   accuracy    to 
about    13  places    and   these   results    were   consistent    with   those    found 
by    the    comparison   program.       In   addition    the   numerical    integrations 
seemed   to   take    somewhat    less    time   than    was    required   by    the  power 
series    method   and   was    not    subject    to  slow   convergence  as    the    ima- 
ginary  part    of    the  argument    became   large. 
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COMPUTER  PROGRAM 


PROGRAM  PHASE       MK  IV       SEPT  1972 
IMPLICIT  REAL*8  (A-H,0-Z) 
INTEGER  RO 
EXTERNAL  CGAR3 , F ACT , DCOS ,DS I N . FQN f FCN 

STCRAGE  FOR  SAMPLE  POINTS  AND  WEIGHTS  FOR  256  AND  512  POINT 
GAUSS-LEGENDRE  INTEGRATIONS. 

COMMON /L25  6/X256( 12 C) , A256(128) 

COMMON/ L512/X5 12 (2  56), A512(256) 

INTEGRATION  INTERVALS  FOR  GAMMA  FUNCTION  INTEGRATIONS 
^/CD/DFC( 15) ,DMC( 15) 

MISCELLANEOUS  CONSTANTS  USED  IN  VARIOUS  SUBROUTINES 
COMMON/ CON ST /PHI ,PH2. FR1, FR2, ASCL, P I ,CA, CT, AA 

PASSES  ALPHA  AND  AND  ETA  TO  THE  INTEGRATION  AND  MATRIX  ROU- 
TINES 

COMMC      *./AL=A,ETA 

PASSES  K  T3  THE  INTEGRATION  ROUTINE 
Cl         /P/WAVN 

PASSES  NUCLEAR  CHARGE 
COMMDN/ATNO/Z 

DUMMY  VARIABLES.   THIS  COMMON  USED  IN  INTEGRATIONS  FOR 
RHO  =  -  1 

COMMON/MCNE/D1 ,D2,D3 

USED  IN  DETERMINING  THE  LIMITS  OF  THE  RANGE  OF  THE  INTEGRAND 
OVER  WHICH  TQ  FIND  THE  ZEROES 

COMMON/I NT/XSM,FF, BO tXLM 

DUMMY  VARIABLES.   THIS  COMMON  USED  IN  INTEGRATIONS  FOR 
RHO  >  -1 

CUMMON/GEZRO/3  4,D5,D6, 110, 120 
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INDICES    FOR    MATRIX    ELEMENTS 

COMMON/ INDEX/M I , NJ , LI , LJ t M I , MJ , IS 

VALUE    OF    EXPONENT    OF    X    IN    BOUND-FREE     INTEGRALS 
COMMON/RHO/RO 

STORAGE    FOR    BOUND-F^EE     INTEGRALS 

DIMENSION    AS (1^52) , 3S { 145  2) , AC (145  2), BC( 14  52) 

STORAGE    FOR    THE    VALJES    3F    THE    EXPONENTS    OF    Rl,     R2    AND    R12     IN 
THE    HYLLERAAS     FUNCTIONS.        LOCATION    OF    THE    REGIONS    DF    ALPHA 
TO    BE    SKIPPED     IN    THE    CALCULATION 

DI  MENS  I  ON    N195  )  ,  L  (  95  )  ,  M  (95  )  ,  ISTPdO  ) 

DATA    INPUTS     FOR    EXA2T     ENERGY     CALCULATIONS 

DIMENSION     X( 5) ,Y(5) 
C  DIMENSION    STATEMENTS    FOR    PHASE    SINGLET,     N=7. 

STORAGE    FOR    BOUND-BOUND    INTEGRALS 
DIMENSION    A(5508) , B(5508) 

STORAGE    FOR    MATRIX    ELEMENTS    AND    EIGENVECTOR    ELEMENTS 
DIMENSION    H(70  ,70)  ,  S  (7  J,7J>),C(70,7J) 

STORAGE    FOR     EIGENVALUES 

DIMENSION    WAV( 70) ,ENERG(70) 

STORAGE  FOR  UPPER  TRIANGULAR  ELEMENTS  OF  H  AND  S  AND  FOR  THE 
INVERSE  OF  THE  TRIANGULAR  FACTOR  OF  S 

DIMENSION  TH(2485) ,TS(2485) 

EQUIVALENCE (H( 1 ) , AS (  1)  ),  (H(2451),BS(  1)  ) 

EQUIVALENCE (S(  1 )  , AC (  1 )  )  , ( S ( 245  1 )  , BC ( 1 ) ) 


DATA  FOR  THE  HYLLERAAS  FUNCTIONS  AND  THE  NUMERICAL  INTEG- 
RATIONS ARE  STORED  IN  A  SEQUENTIAL  DATA  SET 

READ(8,1000)N, L,M 

READ(8,1002)  DPCDMC 

READ(8,10  01)  U512( I  I) , A512(  II ),  11  =  1,25  6) 

READ(8,10  01)(X256(KK),A256(KK) ,KK=1, 12  8) 

REMAINDER  OF  DATA  FIR  PRESENT  CALCULATION  IS  READ  IN  FROM 
PUNCHED  CARDS  SUBMITTED  WITH  THE  MAIN  PROGRAM.   THIS  DATA 
REMAINS  FIXED  FOR  THE  ENTIRE  PROGRAM:   INDEX  FOR  NUMBER  OF 
INTEGRATIONS,  TOTAL  SPIN,  BASIS  SET  SIZE  DETERMINED  BY  THIS 
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INDEX,  NUMBER  OF  ENERGIES  TO  BE  INVESTIGATED  (SET  THIS  EQUAL 
TO  ONE  Ic  THE  ENTIRE  ELASTIC  RANGE  IS  TO  BE  COVERED), 
NUCLEAR  CHARGE 

RE AD { 5, 1021)  NSET, I S t MATRI X , I  NO , Z 

IF( IS.EO.O )  GD  TO  88 

IF    TOTAL    SPIN     IS    ONE     (TRIPLET)     SCREEN    OUT    THOSE     FUNCTIONS 
NOT    NEEDED    FOR    THE    CALCULATION 

K    =    0 

DO     100    KK=1,95 

IFIN(KK) .EQ.L( KK ) )     GO    TO    100 

K    =    K    +     1 

N(K)     =     N(KK) 

L(K)     =     L(KK) 

M(K)     =     M(KK) 

IF(K.  EG.  MATRIX)     KK    =    95 

100    CONTINUE 

SET    VARIOUS    CONSTANTS    WHICH    DETERMINE    NEEDED    INTEGRALS 
88    MAX    =    2 -NSET    +     4 
I  MAX    =    NSET    +     5 
MAXI     =     NSET    +    3 
LIM    =    MATRI X*MATRI X 
I  SIZE    =    MATRIX* (MATRIX    +    l)/2 

FIND    THE    REQUIRED    BOUND-BOUND    INTEGRALS 
CALL    TBLKMAXt AtBJ 

SET    CONSTANTS    NEEDED    FOR    THE    NUMERICAL    INTEGRATIONS 
PI     =    0. 3i41592653589793Dl 
XSM    =    0.5D-2 
FF    =    0.3D-1 
BO    =    0.1oD2 
XLM    =     0.45  5D2 

FIRST    LOOP.       READ    D*TA    FOR    THE    RANGE    OF    ALPHA    AND    ENERGY    TO 
BE     INVESTIGATED. 

DO    130    I W  =  1,1  NO 

RE  AD (5,1 008 )  I ALC,  I  AH  I , I  A I  NT , NUM, FCT , PLD , PHI 

VALUES  OF  ALPHA  TO  BE  SKIPPED 

READ(5,1005)  ISTP 
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USED  TO  LIMIT  SEARCH  FOR  CLOSEST  EIGENVALUES  IN  EXACT  ENERGY 
CALCULATIONS.   IF  ND  EIGENVALUE  IS  FOUND  WITHIN  THIS  EX- 
TENDED RANGE  THE  CALCULATION  IS  STOPPED  FOR  THAT  ENERGY 

PHIEXT  =  PHI  +  O.lD-1 

PLOEXT  =  PLO  -  O.lD-1 

SECOND  LOOP.   COVERS  DESIRED  RANGE  OF  ALPHA  FOR  DESIGNATED 
ENERGY 

1J)2  DO  101  I  A=I  ALO,  IAHI  ,IAINT 

IF  THE  LOW  VALUE  OF  K  DESIRED  IS  ZERO,  THE  ENTIRE  ELASTIC 
scatter;     EGION  WILL  BE  INVESTIGATED.   SKIP  THE  EXACT 
ENERGY  CALCULATION 

IF(PL3.EO. j. 3D3)  GO  TO  50 

READ  THE  VALUES  OF  ALPHA  AND  K  FOR  THE  FIVE  POINTS  ON  THE 
TRAJECTORY  CLOSEST  TO  THE  DESIRED  VALUE  OF  K 

READ (5, 1016) (Y(I)tX(I) .1=1,5) 

XINT  =  PLO  +  0.5D-6 

USING  A  CUBIC  INTERPOLATING  POLYNOMIAL,  FIND  AN  ESTIMATE  OF 
THE  VALUE  OF  ALPHA  REQUIRED  TO  GIVE  THE  DESIRED  VALUE  OF  K 

70  CALL  SPLINH  X,  Y,5,  XINT,YIMT) 

ALFA  =  YII 

GO  TO  51 

CHECK  TC  SEE  IF  A  SEGMENT  OF  THE  ALPHA  RANGE  IS  TO  BE  SKIP- 
PED.  NOT  JSED  WHEN  THE  EXACT  ENERGY  CALCULATION  IS  USED 

50  DO  103  I  =  1,NIW,2 

IF  (IA.EO. ISTP( I ) )  IA=ISTP(I+1) 
103  CONTINUE 

ALFA  =  FCT*DFLOAT(IA) 

EVALUATE  THE  H  AND  S  MATRIX  ELEMENTS 

51  DO  16  J=l. MATRIX 
J  J  =  J* (J  -  1) /2 

DO  10  1=1 , J 
N I  =  N  (  I  ) 
NJ  =  NU) 
LI  =  L( I ) 
LJ  =  L( J) 
M  I  =  M  (  I  ) 

; :  j  =  m  (  j ) 
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CALL  ELEM(HH,SS,A,B,MAX) 
IJ  =  I  +  JJ 
H(  I  ,J)  =  HH 
TH( IJ)  =  HH 
S(  I  ,J)  =  SS 

10  TS( IJ)  =  SS 
16  CONTINUE 

BEFORE  FACTORING  S,  SCALE  IT  SO  THE  LARGEST  ELEMENT  IS  <  10, 
F  I S  THE  SCALc  FACTOR 

CALL  SCUTS,  IS  IZE,  F) 

WRITE( 6, 1006)  F 

CALL  SMPYl TS,F, ISIZE) 

FACTOR  S  INTO  ITS  TRIANGULAR  ELEMENTS  AND  INVERT  THE  UPPER 
TRIANGULAR  FACTOR.   RETURN  IT  IN  TS 

CALL  INVERT  (TS»H,St  CM  ATRIXt  IS  IZE) 

FORM  U**{-T)*H*U**(-1).   RETURN  IT  IN  TH 
CALL  TMPRD(TH»TS,C,MATRIX,F) 

EXPAND  TH  TO  FULL  SIZE  AND  STORE  IN  H.   PUT  TS  IN  UPPER 
TRIANGLE  OF  S  AND  SAVE  FOR  FUTURE  USE 

DO  2  J  J  =  1, MATRIX 

JJ  =  J  *  ( J  -  1  J  /  2 

DO  20  1=1, J 

I  J  =  I  +  JJ 

HT  =  H(  I  ,  J  ) 

H(  I  tJ]  =  TH(I J  ) 

H(  J,I  )  =  TH( IJ  ) 

TH(  I  J)  =  HT 

S( J, I )  =  O.ODD 

20  S(  I , J)  =  TS(IJ) 

TRACE  =  O.ODO 

DO  11  1=1, MATRIX 

11  TRACE  =  TRACE  +  H( I  ,1  J 
WRITE16.1013)  TRACE 

FIND  EIGENVALUES  AND  EIGENVECTORS  OF  TRANSFORMED  H  MATRIX 
CALL  DJACVT(H,MATRIXtl , E NERG , C , MAT RIX ) 
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TRACE    =    O.ODO 

FIND    SUM    OF    EIGENVALUES    AND    COMPARE    WITH    TRACE    OF    H     PREV- 
IOUSLY   FOUND 

DO    17     1=1, MATRIX 

17    TRACE    =    TRACE    +    ENERG(I) 

WRITE(6,1015)    TRACE 

SCALE    EIGENVECTORS    SO    LARGEST    ELEMENT    IS    <    10 
CALL    SCL(C  ,LIM,F) 
CALL    SMPY(C,F,LIMJ 

MULTIPLY    THE    ORTHOGD^AL    EIGENVECTORS    BY    U**(-l)     TO    FIND 
THE    EIGENVECTORS    OF    THE    ORIGINAL    PROBLEM 

CALL    MPRD(S»C» H, MATRIX) 

WRITE(6f 3002J  ALFA 

WRITE(6,3004  )  ENERG 

IK  =  0 

THIRD    LOOP.        EXAMINE    THE     EIGENVALUES    TO    FIND    THOSE     IN    THE 
DES  IRED    ENERGY    r, 

DO    15    I E=lf MATRIX 

WAV( IE)     =    O.ODO 

E    =    ENERG ( IE ) 

ROOT    =0.2D1*E    +    1*1 

IF(ROOT .LE.0.3D0)     ROOT    =    O.ODO 

PINT    =    DSQRTt  ROOT) 
IF(  (PINT.  LE.PLO) .OR.  (PINT. GT  .PHI ) ) GO    TO     21 

INDEX    GIVING    NUMBER    OF    EIGENVECTORS    IN    THE    DESIRED    ENERGY 
RANGE 

IK    =     IK    +    1 

WAV(IK)     =    PINT 

EIGENVALUES    IN    DESIRED    RANGE    RELOCATED    IN    FIRST    ADDRESSES    OF 
ENIGENVALUE    VECTOR 

ENERG (IK)     =    E 

EIGENVECTOR  CORR ESPDNDING  TO  ABOVE  EIGENVALUE  RELOCATED  IN 
APPROPRIATE  COLUMN  OF  EIGENVECTOR  MATRIX 

DO  12  K=l, MATRIX 

12  C(K, IK)  =  C(K, IE) 
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IN  EXACT  ENERGY  CALCULATION  ASSUME  FIRST  VALUE  WITHIN 
DESIRED  LIMITS  THE  THE  ONE  LOOKED  FOR 

IF(PLC.3T.0.0D0)  IE  =  MATRIX 

GO  TO  15 

IF  EXACT  ENERGY  FEATURE  IS  NOT  BEING  USED  SKIP  THIS  SECTION 
21  IF(PL3.EQ. 3.0D3)  GO  TO  15 

IF  EIGENVALUE  BEING  EXAMINED  IS  OUTSIDE  THE  EXTENDED  RANGE 
GO  ON  TO  THE  NEXT  QME 

IF( (PINT. LT.PLOEXT). OR. (PINT. GT.PHIEXT) )GO  TO  15 

IF  EIGENVALUE  IS  OUTSIDE" THE  DESIRED  RANGE  BUT  WITHIN  THE 
EXTENDED  3NE ,  PRINT  THE  VALUE  AND  REPLACE  THE  OUTERMOST 
POINT  ON  THE  SEGMENT  OF  THE  TRAJECTORY  UNDER  EXAMINATION  BY 
THIS  VALUE.   THEN  GO  BACK  AND  RECALCULATE  AN  IMPROVED 
ESTIMATE  3F  THE  REQUIRED  ALPHA 

WRITE (6 , 1014)  IE, PINT 

KNUM  =  0 

IF(PINT.LT.X(1 ) )  KNUM  =  1 

IFCPINT .GT.XC5  J  J  KNUM  =  5 

IF( (KNUM.EO.l) .OR. (KNUM.EQ.5) )  GO  TO  65 

DO  60  I =  2, 5 

IF((PINT.LT.X(  I)  ). AND. (PINT. GT.X(  1-1))  J  KNUM  =  I 

6J  CONTINUE 

IF{ (KNUM.LT.l ) .OR. (KNUM.GT.5) )  GO  TO  15 

IF(  ( PINT-X(  1))  .GE. ( XI 5)-PINT)  )  GO  TO  61 

DC  62  1=1,4 

IL  =  5  -  I 

IF(  IL.GE.KNUM)  XUL+1)  =  X(IL) 

IF(  IL.GE.KNUM)  YUL+1)  =  Y(IL) 

62  CONTINUE 
X(KNUM)  =  PINT 
Y(KNUM)  =  ALFA 
GO  TO  70 

61  KNUM  =  KNUM  -  1 
DO  63  1=2,5 

IF(I.LE.KNUM)  X { I - 1 )  =  X( I ) 
I  F (  I .LE.KNUM)  Y(  1-1)  =  Y(  I  ) 

63  CONTINUE 

65  X (KNUM J  =  PINT 
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Y(KNUM)  =  ALFA' 
GO  TO  70 
15  CONTINUE 

IF  NO  EIGENVALUES  WERE  FOUND  IN  THE  DESIRED  RANGE,  GO  TO  THE 
NEXT  VALUE  OF  ALPHA 

IF(IK. EO.O)  GO  TO  101 

WRITE(6,3005)  ALFA 

WRITE (6,3004) (WAV ( LP) ,LP=1 , IK) 

FOURTH  LOOP.   BEGIN  THE  PHASE  SHIFT  CALCULATION 
DO  113  NN=1, IK 
E  =  ENERGCNN) 
WAVN  =  WAV(NN) 

FIRST  FIND  THE  COULOMB  PHASE  SHIFT,  ETA 
AA  =  (Z  -  0.1DD/WAVN 

IF  Z  =  1  SKIP  THIS  PART  AND  GO  ON  TO  THE  BOUND-FREE  INTEG- 
RATIONS 

IF(AA. EO.O. ODD )  GO  TO  301 

ASO  =  AA*AA 

CA  =  ASQ    -  0.11D2 

CT  =  0.6D1*(ASQ    -  0.1DD/AA 

CALL  INT256(C3ARG,GR) 

CA  =  CT 

CT=  0.11D2  -ASO 

CALL    INT256(C3ARG,GI ) 

GR    =    -GR 

ETA    =    DATAN2(GI ,GR) 

GO  TO  302 
301    ETA    =    O.ODO 

WRITE (6,1323)     MATR I  X , M ATR I  X , AL FA , E,WAVN 

GO  TO  303 
332    WRITE(6,1004)     MA TRI X , M AT RI X , ALF A, E , WAVN , ET A, Z 

FIND    THE    REQUIRED    BOUND-FREE     INTEGRALS 
303    CALL    TBL2( IMAX,MAXI , AS , AC , BS , BC ) 
STERM    =    O.ODO 
CTERM    =    O.ODO 
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FORM    THE    BOUND-FREE    MATRIX    ELEMENTS 
DO    300    JM=1, MATRIX 

CC    =    C(JM,NN) 
N  J    =    N  (  J  M ) 
LJ    =    L(JM) 
M  J    =    M  (  J  M ) 

CALL     PHASE (AS,  AC , BS , BC , MAX  J .  , HS , HC , E J 
HSC    =    HS*CC 
STERM    =    STERM    +    HSC 
HCC    =    HC*CC 
CTERM    =    CTERM    +    HCC 
300    CONTINUE 

EVALUATE    THE    PHASE    SHIFT 

TDLTA    =    -STERM/CTERM 
TDLTAK     =    TDLTA/WAVN 
DLTA    =    DATAN(TDLTA) 
IF ( IS.EQ.O)     G3    TO    91 
WRITE (6,  1011  ) 
GO    TO    92 

91  WRITE16.1012 J 

92  WRITEC 6,  1007)  ALFA , WAVN , TDLTA , TDLTAK, D LT A 
WRITE (7, 11 03 J  ALFA, WAVN, T DLT A, T DLT AK , DLT A 
IF(DLTA.GT.O.ODO)     GO    TO    113 

IF    THE    PHASE     SHIFT    IS    NEGATIVE,     ARBITRARILY    ADD     PI     TO    IT    AND 
PRINT     AN    ADDITIONAL    DATA    CARD    WITH    THE    NEW    PHASE     SHIFT 

DLTAP    =    DLTA    +    PI 

WRITE(7,1101 )     ALFA, WAVN, DLTAP 
113    CONTINUE 
131    CONTINUE 
130    CONTINUE 

STOP 

1000  F0RMAT148I 1,/,47I1  ) 

1001  FORMAT (2D30. 16 ) 

1032  F0RMAT(8D9.3,/  ,7D9.3) 

1003  FORMAT (•-• ,5X, ' SOLUTIONS  FOR  ALPHA  = • , Fl 1. 6t /t •- • J 
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1004  FORMAT (  '  1'  t /, '-•  , lOXt  '  ELECTRON  -  HELIUM  +  SCATTERING 
1PHASE  SHIFT  FOR'  ,13, •  X  ',12,'  HAMILTONIAN  MATRIX.',/, 
2'-', 34X,  "ALPHA  ='  ,F10.  7,/, « 0'  ,24X, 'TOTAL  ENERGY,  E  =  ', 
3D24. 16,/, «0'  ,25X, • WAVE  NUMBER,  K  =  •  , D24. 16  ,  /  ,  • 0'  , 1 5X , ' 
4CCULCMB  PHASE  SHIFT,  ETA  =  '  , D24 . 16 ,/ , ' 3 « , 23X ,  •  "NUCL EAR 
5CHARGE",     Z     =',F9.5) 

1005  FORMAT(10I8) 

1006  FORMAT(  • 1'  ,  10X,  'F  =  « , D 10. 3 , /, • - ■ ) 


1037  FORMAT(  •  0 


:X, 'ALPHA  =•  ,F9.6,4X, 'K  =•  , D24 . 16 , 5X, ' 


1TANGENT (PHASE  SHIFT)  = ' , D24. 16 , / , » 0 { , 65X , ' TAN( DE LTA ) /K 
2='  ,D24.16,/,' 3'  ,66X, 'PHASE  SHIFT  =',D24.16) 

1008  FORMAT     ( 8X , 21 8 , 2  I  3, F 10. 3 , 5X , 2F 1 5. 9) 

1011  FORMAT { ' 3' ,/,'-' ,10X, « TRIPLET  PHASE  SHIFTS.') 

1012  FORMAT! '0'  ,/,'-'  ,10X, 'SINGLET  PHASE  SHIFTS.') 

1013  FORMAT( 10X, 'TRACE  OF  f U**(-T ) ) *H*(U**(-1 ) )  =',D24.16) 

1014  FORM  AT C '-«  ,10X,'  K( •  ,12  ,  '  )  =  »  ,024 . 16 ,/,'-•  ) 

1015  F0RMAT(22X, ' SUM    OF    EIGENVALUES    =',D24.16,/) 
10  16  F0RMAT{F11.7,F12.7) 

1C17    FORMAT {'-'  ,10K  ,  15,  '    CARDS    PUNCHED    FOR    THIS    JOB.') 

1020  FORMAT ( ' 1' ,/,•-« ,10X, '     ELECTRON    -    HYDROGEN     SCATTERING 

1  PHASE  SHIFT  FOR ',13,'  X  ',12,'  HAMILTONIAN  MATRIX.',/, 
2'-'  ,34X,  'ALPHA  = '  , F 9. 6 , / , • 0'  , 24X , ' TOT AL  ENERGY,  E  =  ', 
3D23.16,/ , '0' ,25X, 'WAVE    NUMBER,    K    =',D23.16) 

1021  FORMAT! 212,214,08.3) 

1100  FORMAT (Fll .7.F12 .7,3X, 3D18.10) 

1101  FORM AT (F 11. 7, F 1 2. 7 , 39X , D 18 . 10) 

3032    FORMAT ('-•  ,5X, • EIGENVALUES    OF    H    FOR    ALPHA    =     «,F12.7,/) 

3004  FORMAT (4X.5D25.16) 

3005  FORMAT  (  '-•  , /,5X,  'THE  VALUES  OF  K  FOR  ALPHA  =',F12.7,' 
1  ARE: ' ,/,' 0' ) 

END 
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SUBROUTINE    T  B.  ].  I  MAX,  A  ,  B  ) 

TBL1    GENERATES    THE    TABLE    OF    BOUND-BOUND    INTEGRALS 
IMPLICIT    REALMS     (A-H,0-Z) 
DIMENSION    A( 1) ,B(1) 

THE     INTEGRALS    ARE    STORED     IN    A    ONE    DIMENSIONAL    ARRAY    CORRES- 
PONDING   TO    THE    VARIOUS    VALUES    OF    NU,     LAMDA    AND    MU 

MAX2    =    MAX*MAX 

LIM    =     (MAX    -    1 J*MAX2 

INITIALIZE    THE    A    AND    B    ARRAYS 
DO    1    1=1, LIM 
A(  I)     =    . 3DD 

i   em   =   .odo 

NMAX    =    MAX    -    1 

EVALUATE    THE     INTEGRALS    FOR    MU    =    2 
DO    2    NN=ltNMAX 
DO    2    LL=2,MAX 

IFUNN+LL.LT.5)  .OR.  CNN+LL.GT  .MAX+2  )  J    GO    TO    2 
IN    =    NN    +    MAX*(LL    -    1)     +    2*MAX2 
CALL    SUM(NN,LL ,AA) 
A(IfJ)     =    AA 

2  CONTINUE 

USE    THE    RECURSION    RELATIONS    TO    FIND    THE     INTEGRALS    FOR    MU    = 
DO    3    NN=2,MAX 
DO    3    LL=2,MAX 
IFUNN+LL.GT.MAX+3)  . OR. (NN+LL .LT .6  }  )    GO    TO    3 

IN    =    NN    +    MAX*(LL    -     1)     +    MAX 2 

INP    =     IN    +    MAX  2    -    1 

A  (  I  N  )     =     A  (  I N  P  ) 

IF((LL.LT.3).0R.(NN.LT.3).0R.(NN+LL.LT.6).0R.(LL.GE. 
1MAX) .OR. (NN.GE.MAX)  }    GO    TO    3 

INP    =    IN    +     MAX* (MAX    +    1)     -    2 

B( IN)     =    A( INP)/ .3D1 

3  CONTINUE 

MMAX    =     MAX    -    I 
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NLMAX    =    MAX    -    2 

USE    THE    RECURSION    RELATIONS    TO    FIND    THE     INTEGRALS     FOR    MU    >    2 

DO    5    MM =4, MM AX 

DO    5    NN=2t NLMAX 

DO    5    LL=2, NLMAX 

IFC (NN+LL.LT.5).OR.(NN+LL+MM.GT.MAX+5)  J    GO    TO    5 

IN    =    NN    +    MAXMLL    -    1)     +     (MM    -    1)*MAX2 

INN    =     IN    -    2*(MAX2    -    1) 

INL    =     IN    -    2*MAX*(MAX    -     1) 

INNL    =    IN    -    MAX* (2*MAX    -    1)     +    1 

A(IN)     =    AC  INN]     +    A(INL)    -     .2D1#B(INNL) 

IF(  ( MM. 3  T.MAX- 3) .0R.(NN+LL.LT.6).CR.(LL.LE.2).0R.(NN. 
1LE.2) .OR.ILL.GT .MAX-3) .CR. ( NN.GT.MAX-3) )     GO    TO    5 

B(IN)     =     (B(INN)     +    B(INL)     -    . 2D1* A( INNL) ) *DFLOAT (MM-3 )/ 
1DFLOAT (MM+1) 

5    CONTINUE 

RETURN 

END 


SUBRCJT  INE    SUVMNU,LMDA,CUT  ) 

SUM    EVALUATES     A( NU ,L AMDA , 2)     FOR    GIVEN    NUt     LAMDA 
IMPLICIT    REAL*8    (A-H,0-ZJ 
EXTERNAL    FACT 
A    =     O.ODO 
NN    =    NU    -    1   • 
DO    10     I-ltLMDA 
IN2    =    LMDA    -    I 
INI     =    NN    +    IN2 
10    A    =    A    +    FACT(IN1)/(FACT(IN2)*(0.2D1**(IN1    +    1))) 
LL    =    LMDA    -    1 

CUT    =    FACT(LL)*(FACT(NNJ    -    A) 
RETURN 
END 
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DOUBLE    PRECISION    FUNCTION    FACT(N) 

FACT     FORMS    N    FACTORIAL 

IMPLICIT    REAL* 8     (A-H,0-Z) 
FACT    =     .1D1 
IF(N.EC.O)     GO    TO    2  0 
DO    ID    1=1, N 
FI    =     I 
10    FACT    =    FACT*FI 
20    RETURN 
END 
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SUBROUTINE  SPL  I  Nl  (  X  ,Y  ,  M  ,X  I  NT  ,  Y  I  NT  ) 

SPLI.N  PROVIDES  INTERPOLATED  VALUE  OF  THE  ORDINATE  USING 
"CUBIC  SPLINE  FITTING" 
USAGE: 

FIRST  CALL  TO  SUBROUTINE: 

CALL  SPLIN1 (X, Y,M,XINT,YINT) 
SJBSEG'JENT  CALLS  USING  THE  SAME  DATA 

CALL  SPLINN(X, Y  ,  M,  X  INT  ,  Y  INT  ) 
DESCRIPTION  OF  PARAMETERS 

X:  MONCTONICALLY  INCREASING  ABSCISSA  ARRAY 
Y:  ONE-FOR-ONE  CORRESPONDING  ORDINATE  ARRAY 
M:  NUMBER  OF  X  AND  Y  VALUES  SUPPLIED  <  OR  = 
300.  (MODIFIED  IN  THE  PRESENT  CONTEXT  TO 
NO  MORE  THAN  FIVE) 
XINT:  VALUE  OF  ABSCISSA  FOR  WHICH  CORRESPOND- 
ING ORDINATE  IS  TO  BE  INTERPOLATED  (OR 
EXTRAPOLATED) 
YINT:  INTERPOLATED  (OR  EXTRAPOLATED)  ORDINATE 
VALUE 
MATHEMATICAL  METHI   : 

UPON  FIRST  ENTRY  TO  SPLIN,  A  CALL  TO  SPLICO  IS 
MADE  TO  DETERMINE  THE  COEFFICIENTS  TO  BE  USED  IN 
PERFORMING  THE  INTERPOLATIONS.   SEARCH  FOR  BRACKET- 
ING ABSCISSA  VALUES  IS  ALWAYS  MADE  FROM  THE  REFER- 
ENCE LAST  JSED  IN  INTERPOLATING. 
REFERENCE: 

PENNINGTON,  R.H.,  "INTRODUCTORY  COMPUTER  METHODS 
AND  NUMERI2AL  ANALYSIS",   MACMILLAN,  NEW  YORK,  1965 

IMPLICIT  REAL<=8  (A-H),REAL*8  (O-Z) 

D I  MENS  ION  X(M) ,Y(M) ,C( 4,6) 

CALL    SPLICOtX, Y,M,C) 

K=l 

EUr  Y  SPLINNU,  Y,M,  XINT,  YINT) 

3  IF (XINT -X( 1)1  70,1,2 

70  K=l 

GO  T3  7 

1  YINT=Y(1) 
RETURN 

2  IF(XINT-X(K+1) )6,4,5 

4  YINT=Y(K+1) 
RETURN 

5  K=K+1 

IF(M-K)  71,71,3 

71  K=M-1 
GO  TO  7 

6  IF( XINT-X( K)  ) 13  ,12,11 
12  YINT=Y(K) 

RETURN 
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13  K=K-1 

GO  TO  6 

7  PRINT  101, XI  NT 

101  FORMAT! 8H0XI\T  =  E18.9,32H,  OUT  OF  RANGE  FOR 
1  INTERPCLAT  ION) 

11  YINT=(X(K+1)-XINT)*(C( 1,K)*(X( K+l ) -XI  NT ) ** 2+C (3 , K)  ) 

YINT=YINT+ (XIMT-X(K) J* (C (2 ,K )* (X INT-XCK ) )**2+C(4,K) ) 

RETURN 

END 


SUBROUTINE  SPL  I  CO ( X , Y , M, C ) 

IMPLICIT    REAL*8     (A-H),REAL*8     (0-Z) 

DIMENSION    X(M)  ,  Y  (  M  )  ,  C  (  4  ,  6  )  ,  D  (  6  )  ,  P  (  6  )  ,  E  (  6  )  ,  A  (  6,  3  )  ,  B  (  6  )  r 
1Z(6) 

MM=M-1 

DO    2    K=l, 

D(K) =  Xl K+l )-X( K) 

P(K)=D(K)/6. 

2  E(KJ=( YtK+l)-Y(K) ) /D(K) 
DC    3    K=2,MM 

3  B(K)=E(K)-E(K- 1) 

At  1,2) =-l.-D( 1) /D(2) 

A(l  ,3 )=Dtl  }/D( 2  ) 

A(2,3)=P(2)-P( l)*A( 1,3) 

A(2,2)=2.*(P(l)+P(2))-Ptl)*A(l,2) 

A(2.3)=A(2,3)/A(2,  2) 

B(2)=B(2)/A(  2,2) 

DO    4    K=3,MM 

A(Kt2)  =  2.*(PU-l)  +  P(K)  )-P(K-l)*A(K-l,3) 

CtK) =B(K)-P [K-1)*B(K-1 ) 

A(K,3)=P(K )/A(  K,2) 

4  B(K)=b(K)/A(K,  2) 
0=D(M-2) /DtM-1 ) 
A(M,  1)=1.+G+A( M-2,3) 

Al M,2)  =-0-A(M,  1)*A(M-1 ,3) 
B(M)=B(M-2 )-A(M, 1)*B{M-1 J 
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Z(M)=B(M)/A(M,  2) 

MN=M-2 

DO  6  I=1,MN 

K  =  M-I 

6  Z(K)=B(K)-ACK,  3)*Z(K+1) 

Z(l) =-A(l,2)*Z 12)-A(1,3)*Z(3) 

DO  7  K=lfMM 

0=1. /( 6.*D(K) ) 

C(1«K)=Z(K)*Q 

C(2tK)=Z(K+l)*Q 

C(3»K)=YCK)/D(K)-Z(K)*P(K) 

7  CC4tKJ=Y(K+l)/D(K)-Z(K+l)*P(KJ 
RETURN 

END 
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SUBROUTINE    ELE^U  HH, SS, A, B, MAX) 

ELEM    GENERATES    THE    MATRIX    ELEMENTS    <CHI ( I  J  | H I  CHI ( J) >    AND 
<CHI  (I )  I  CHI ( J)>    FCR    H    AND    S,     RESPECTIVELY 

IMPLICIT    REAL*8     (A-H.O-Z) 

INTEGER    PWR 

COMMON/TRM/ALFA, ETA 

COMMON/ATNO/Z 

COMMON/ INDEX/NI , NJ,LI,LJ ,MI,MJ,  IS 

DIMENSION    A( 1)  , 3(1 ) ,NU( 52) ,MU( 52) tLMDA( 5  2)  ,COEF( 52) 

NIJ    =    NI    +     NJ 

LIJ    -     LI    +    LJ 

NLIJ    =    NI     +    LJ 

LNIJ    =    LI     +    NJ 

MIJ    =    MI    +    MJ 

FN1     =    DFLOAT(^JJ     +    1) 

FL1    =    DFLOAKLJ    +    1) 

FN    =    DFLOAT (NJ ) 

FL    =    DFLOAT (LJ  ) 

FM    =    DFLOAKMJ  ) 

MAX2    =    MAX*MAX 

SET    THE    VALUES    OF    NJ  ,     LAMDA,     MU    AND    THE    COEFFICIENTS    FOR 
EACH    OF    THE    52    TERMS    OF    H(I,J)     AND    THE    4    TERMS    OF     S(ItJ) 

NU(     1)    =  NIJ    ♦    2 

NU(  3)  =  NU(1) 

NU(  4)  =  NU(  1) 

NU(  6)  =  NU( 1) 

NU(  8)  =  NU(1) 

NU(  9)  =  NU( 1) 

NU(37)  =  NU(1) 

LMDA(IO)  =  NU(  1) 

LMDA(ll)  =  NIM  1) 

LMDA(14)  =  NU(  1) 

LKDA115)  =  NUJ  1) 

LMDA(16)  =  NU(  1  ) 

LMDA118  )  =  NU(  1  ) 
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LMDA(42) 
NU(  2)  = 
LMDAl 12) 
NU(39)  = 
LMDA(44) 
NU(  5)  = 
NU138)  = 
NU(40)  = 
LMDAC131 
LMDA141 ) 
LMDA(43) 
NU(  7)  = 
LMDA( 17  ) 
NU(10)  = 
NUC12J  = 
NU(13)  = 
NU(15)  = 
NU(17)  = 
N  U  (  1  8 )  = 
NU(41)  = 
LMDA(  1) 
LMDA(  2) 
LMDAl  5) 
LMDA(  6) 
LMDAl  7) 
LMDA(  9) 
LM0A{38) 
NU(ll)  = 
LMDA(  3) 
NUC43J  = 
LMDA140 J 
NU(14)  = 
LMDA(  4) 
NU(42)  = 
NU(44)  = 
LMDAC37) 


=  NU(  1) 
NIJ  +  1 
=  NU(  2) 
NU(2) 
=  NJl  2) 
NIJ  *  3 
NU15) 
N  U  ( 5 1 
=  NU(  5) 
=  NU( 5) 
=  NU(  5) 
NIJ 

=  NU<  7) 
LI  J  v  2 
NU ( 10 1 
N  U ( 10  J 
NUC1Q ) 
NU(IO) 
NU(  10) 
NU ( 10  J 
=  \'U(  10) 
=  NU( 10) 
=  Nll(  10) 
=  NU( 10) 
=  NU(IO) 
=  NU( 10) 
=  NU( 10) 
LIJ  +  1 
=  NUl  11  ) 
NUI11  ) 
=  NU{  11  ) 
LI  J  *  3 
=  NU( 14) 
NU( 14) 
NUC14  ) 
=  NUl 14) 
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LMDAC39) 
NU(16)     = 
LMDA(     8) 
NUU9J     = 
NU(21)     = 
NU(22)     = 
NU(24J     = 
NU(26)     = 
NU(27)     = 
NU(45)     = 
LMDA(28) 
LMDA(29  ) 
LMDA(32) 
LMDA(33) 
LMDA(34) 
LMDA(36) 
LMDA(50 J 
NU(20)     = 
NU(47)      = 
LMDAC30  ) 
LMDAC52) 
NU(23)     = 
NU146)     = 
NU(48)     = 
LNDAC31 ) 
LMDA(49J 
LMDA(51 ) 
NUC25I     = 
LMDA(35 J 
NU(28)     = 
NU(30)     = 
NUC31)     = 
NUC33J     = 
NU(35)     = 
NU(36)     = 
NU(49)     = 


=    NU(14) 
LIJ 

=    NU( 16) 
NLIJ     +    2 
NU( 15) 
NU(19) 
NU(  19  ) 
NU( 19) 
NU(19  ) 
NU(  19) 
=    NU( 19) 
=    NU(19) 
=     NU( 19) 
=    NU(19) 
=    NU(  19  J 
=    NU( 19) 
=     KU(19) 
NLIJ     +     1 
NUC23  ) 
=    NU( 20) 
=    NIK20) 
NLIJ     +    3 
NU(23) 
NU(23  ) 
=    NJ( 23) 
=    NU(23) 
=    NU(23) 
NLIJ 

=  NJ(25) 
LNIJ  +  2 
NU(23 ) 
NU(2B  ) 
NU(28) 
NU(23  ) 
NU(28  ) 
NU( 23 ) 
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LMDA(19  ) 
LMDA(20) 
LMDA{23) 
LMDA(24) 
LMDA(25) 
LMDA(27 ) 
LMDA146) 
NU(29)  = 
NU(51)  = 
LMDA(2  1) 
LMDAK8) 
NU(32)  = 
LMDA122) 
NU(50)  = 
NU(52)  = 
LMDA(45 ) 
LMDAI 47) 
NU(34)  = 
LMDAC26) 
MU(     1)     = 

;    2)    = 

.     3)  = 

MU(     7)  = 

MU  (     8  )  = 

MU(1D)  = 

MU(ll)  = 

MJ(12)  = 

MU(16)  = 

MU(17)  = 

MUC191  = 

MU(20)  = 

MU(21)  = 

MU(2  5)  = 

MU(26)  = 

MU(28)  = 

MUC29J  = 


=  NU128) 
=  NU( 28) 
=  NU  (  2  8  ) 
=  NUC28) 
=  NU(28) 
=  NU( 28  ) 
=  NU( 28) 
LNIJ  +  1 
NU(29  ) 
=  NUC29) 
=  MU( 29  ) 
LNIJ  +  3 
=  NUC32J 
NU(32 ) 
NU(32) 
=  NU(32  ) 
=    NU(32) 

LNIJ 

MU( 34) 
MU    +     2 

MU  (  1  ) 
MU(  1) 

MUll) 

MU  (  1 ) 

MU(  1) 

MU  (  1  ) 

MU(  1) 

MU(1) 

MU(  1) 

MU(  1) 

MU  (  1  ) 

MU(  1) 

MU(1) 

MU(  1) 

MU(  1) 

MU  (  1  ) 
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MU(30)  =    MUll) 

MU{34)  =    MU(1) 

MU(35)  =    MU(1) 

MU(     4)  =    MI  J 

MU(     5)  =    MU(4) 

HUC     6)  =    MU(4) 

MU(13)  =    MU(4) 

MU(14)  =    MU (4) 

HUC 15]  =    MU( 4) 

MU( 22  J  =    MU(4) 

MUC23]  =    MU(4) 

MU(24)  =    MU(4) 

MU(31)  =    MU(4) 

MU(32)  =    MUC4) 

MU(33)  =    MU(4) 

MU(37)  =    MU(4) 

MU(38)  =    MU(4) 

MU(39J  =    MU(4) 

MU(4C)  =    MUC4J 

MU ( 4  1  )  =    MU ( 4 ) 

MU(42)  =    MU(4) 

MU143J  =    MU(4) 

MU(44)  =    MU(4) 

MU(45)  =    MU(4) 

MU(4C)  =     MU(4) 

MU(47)  =    MU(4) 

MU(48)  =    MU(4) 

MU(49)  =    MU(4) 

MU(50)  =    MU(4) 

MU(51)  =    MU(4) 

MU(52)  =    MU(4) 

MU(     9)  =    MI  J    +     1 

MU(18)  =    MU(9) 

MU{2  7)  =    MU(9) 

MU(3  6)  =    MU(9) 

COEF(     1)     =    -0. 25D0*ALFA*ALFA 
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COEF( 10)  =  COEF( 1) 

CCEF(19)  =  COEF(l) 

CCEF(28)  =  COEF(l) 

COEF{     2)  =  0.5DD*ALFA*FN1    -    Z 

COEF(12)  =  COEF(2) 

COEF(21J  =  COEF(2) 

COEF(29)  =  COEF(2) 

COEF(     3)  =  0.5D0*ALFA*FL1    -    Z 

COEF(ll)  =  COEF(3) 

COEF(20)  =  C0EF(3) 

COEF(30)  =  COEF(3) 

COEFt     4)  =  0.5DO*ALFA*FM 

COEF(     5)  =  C0EF{4) 

COEFQ3)  =  C0EF(4) 

C0EF(14)  =  C0EF(4) 

COEF(22)  =  C0EF(4) 

CCEF(23)  =  C0EF(4) 

C0EF(31)  =  C0EF(4) 

CCEF(32)  =  C0EF(4) 

CCEF137)  =  -C0EF(4) 

CCEF(38)  =  -COEF (4) 

C0EF(41 )  =  -C0EF(4) 

COEF(42)  =  -C0EF(4) 

COEF (45)  =  -COEF (4) 

C0EF(46)  =  -C0EF(4) 

C0EF(49)  =  -COEF (4) 

COEF(50)  =  -C0EF(4) 

COEFl     6)  =  -FM*(FM    +    FN    +    FL1) 

COEF(15)  =  C0EF(6) 

C0EF(24)  =  C0EF(6)     . 

COEF{33)  =  C0EF(6) 

COEF(    7)  =  -0.5D0*FN*FN1 

COEF(17)  =  C0EF(7) 

C0EF(26)  =  C0EF(7) 

C0EF(34)  =  COEF(7) 

COEFt     8)  =  -3. 5D0*FL*FL1 
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C0EF(16)    =     C0EF(8) 
CCEF(25)     =    C0EF(8) 
C0EF(35)     =    C0EF(8) 
COEFt     9)    =    0.1D1 
COEF(18)     =    COEF(9) 
C0EF(27)    =    C0EF(9) 
C0EF(36)     =    COEF(9) 
C0EFO9)     =    FM*FN 
CCEF(44)     =    C0EFO9) 
C0EF{48)     =    C0EFO9) 
COEF(51)    =    C0EF(39) 
C0EF(40)     =    FM*FL 
C0EF(43)     =    C0EF140) 
CGEF(47)     =    CCEF140) 
COEF(52)     =    C0EF(43) 
MAX2    =    MAX*MAX 
riH  =    O.ODO 

IF(  IS.EO.O)     GO    TO    10 

FOP    TRIPLET    SCATTERING    THE    SIGNS    OF     THE    LAST    26    TERMS    IN 
H( I ,J)     MUST    BE    CHANGED 

DO    1     1=19,36 

1  COEF( I )     =    -CGEFl I ) 
DO    2     1=45,52 

2  COEF(  I  )     =    -CCEFt  I  ) 

THE    VALUES    OF    NU ,     LAMDA    AND    MU    FOUND    ABOVE    MUST    BE    CONVERTED 
TO    A    SINGLE     INDEX    CORRESPONDING    TD    THE     STORAGE    MODE    OF    THE 
INTEGRALS 

10     DO    6    K=l,36 

IN    =    NU(K)     +    1     +    MAX-LMDA(K)     +    MU(K)*MAX2 

PWR    =    NU(K)     +    LMDA(K)     +    MU(K) 

6  HH  =    HH  +-    COEF(K)*A(IN)/ALFA**PWR 
DO    7    K=37,52 

IN    =    NU(K)     +    1     +     MAX*LMDA(K)     +    MU(K)*MAX2 
PWR    =    NU(K)     +    LMDA(K)     +    MU(K) 

7  HH  =     HH  +    COEF(K)*B( IN)/ALFA**PWR 

AFTER     FORMATION    OF    H(I,J)     THE    SAME    PROCESS    MUST    BE    REPEATED 
FOR    S( I ,J) 
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SS    =    3.0D0 
DO    9    K=l ,10,9 

IN    =    NU(K)     +    1     +    MAX*LMDA(K)     +    MU(K)*MAX2 
PWR    =    NUCK)     +    LMDA(K)     +    MUCK) 
9     SS    =    SS    +    AC  TJ)/ ALFA** PWR 
IF     (IS. EQ.O)     30    TO    II 
DO    13     K=19,28,9 

IN    =    NU(K)     +    1     +    MAX*LMDACK)     +    MU(K)*MAX2 
PWR    =    NUCK)     +     LMDA(K)     +    MUCK) 
13    SS    =    SS    -    AC  INI  )/ALFA**PWR 
GO    TO    12 

11  DO    8    K=19,28,9 

IN    =    NUCK)     +    1     +    MAX-LMDA(K)     +    MU(K)*MAX2 
PWR    =    NUCK)     +    LMDACK)     +    MUCK) 
8    SS    =    SS    +    A(IN)/ALFA**PWR 

12  RETURN 
END 
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SUBROUTINE  SCL(S,N,F) 

SCL  FINDS  THE  ORDER  OF  MAGNITUDE  OF  THE  LARGEST  (ABSOLUTE 
VALUE)  ELEMENT  IN  A  GIVEN  MATRIX 

IMPLICIT  REALMS  (A-H»0-Z) 

DIMENSION  S(l) 

TCP  =  DABS(S(1  )  ) 

DO  100  I=1,N 

T  =  CABS(S (  I  )  ) 

IF(T.GT.TOP)  TOP  =  T 

100  CONTINUE 

FF, =  DLOGIO(TDP) 

IPWR  =  IDINT(FF) 

F  =  0.1D2**(-IPWR) 

RETURN 

END 


SUBROUTINE  SM3Y(A,C,N) 

SMPY  MULTIPLIES  A  MATRIX  BY  A  SCALAR,  RETURNING  THE  SCALED 
MATRIX  IN  THE  SAME  LOCATION  AS  THE  ORIGINAL  MATRIX 

LICIT  REAL*8  IA-H.O-ZJ 

DIMENSION  A(l) 

DO  1  1=1, N 

1  All)  =  MI  )*C 

RETURN 

END 
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SUBROUTINE  IN\/ERT{TS,H,S,C,N,ISIZE) 

INVERT  TAKES  THE  SYMMETRIC,  POSITIVE  DEFINITE  MATRIX  S 
(STORED  IN  THE  COMPRESSED  MODE  IN  TS)  AND,  USING  DSINV,  FAC- 
TORS IT  INTO  ITS  UPPER  TRIANGULAR  FACTOR  AND  INVERTS  THIS 
FACTOR.   THEN  USING  THIS  RESULT  AS  A  FIRST  APPROXIMATION  THE 
INVERSE  IS  IMPROVED  BY  AN  ITERATION  METHOD 

IMPLICIT  REAL*8  (A-H,0-Z) 

REAL*4     EPS 

DIMENSION    TS(1  )  ,H( 1)  ,S(  1)  ,C{  1) 

KIT    =    0 

EPS    =    0.1E-7 

CALL    DSINV(TS,C ,N,EPS, IER) 

IFI IER.EQ.-l )     GO    TO    16 

IS  IG    =     1 

U**(-1J  (IN  TS)  AND  U  (IN  C)  ARE  MULTIPLIED  TO  FORM  A  TRIAL 
IDENTITY  rtHlCH  IS  CHECKED  FOR  OFF-DIAGONAL  ELEMENTS  WHICH 
MAY  BE  >  10**(-20)i  AND  STORED  IN  THE  LOWER  TRIANGLE  OF  H 

DO  3  K  =  2,N 

KK  =  K  -  1 

KM  =  K^:(K  -  1)  /2 

DO  2  1=1, KK 

I  K  =  K  +  N*  (  I  -  1  ) 

CC  =  O.ODO 

DO  1  J  =  I  ,K 

IJ  =  J*( J  -  1 ) /2  +  I 

JK  =  KM  +  J 

1  CC  =  CC  +  C( IJ )*TS (JK) 
IF(DABS(CC ) .GT .0.1D-20)  ISIG  =  0 

2  H( IK)  =  -CC 

3  CONTINUE 

IF  THE  INVERSE  IS  ALREADY  GOOD  ENOUGH  NO  IMPROVEMENT  IS 
NEEDED  AND  RETURN  TO  THE  MAIN  PROGRAM  CAN  BE  MADE 

IF( ISIG.EO.l)  GO  TO  15 

14  KIT  =  KIT  +  1 

THE  CORRECTION  TERM  IS  FORMED  BY  MULTIPLYING  U**(-l)  AND  THE 
TRIAL  INVERSE 

DC  6  K=2,N 


130 


KK  =  K  -  1 

KM  =  K*(K  -  1)  /2 

DO  5  I=1,KK 

IN  »  N*(I  -  II 

IK  =  K  +  IN 

11=1+1 

CC  =  0.000 

CO  4  J=II,K 

IJ  =  J  +  I  N 

JK  =  KM  +  J 

4  CC  =  CC  +  H(IJ  J*TS ( JK) 

5  S( IK)  =  CC 

6  CONTINUE 

THE  IMPROVED  INVERSE  IS  FORMED  BY  ADDING  THE  CORRECTION  TERM 
TO  U**C-11 

DO  8  K  =  2,N 

KK  =  K  -  1 

KK  =  K*(K  -  1)  /2 

DO  7  I =1 , KK 

IK  =  KM  +  I 

KI  =  K  +  N*  (  I  -  1) 

7  TS( I<)  =  TS(  K )  +  S(KI  ) 

8  CONTINUE 
IS  IG  =  1 

THE  NEW  TRIAL  INVERSE  IS  FORMED  AND  TESTED  FOR  ANY  ELEMENTS 
>  10~« (-20).   THIS  VALUE  IS  PUT  IN  THE  LOWER  TRIANGLE  OF  S 

DO  11  K=2,N 

KK  =  K  -  1 

DO  10  1=1,  KK 

IN  =  N*(I  -  1) 

IK  =  K  +  I N 

CC  =  O.ODO 

11=1+1 

I  F (  I  I  .GT .KK)  30  TO  10 

DO  9  J=II ,KK 

IJ  =  J  +  IN 
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JK    =     K     +    N*(  J     -     1) 
9    CC    =    CC    +    H( IJ )*H(JK) 

IF(DA3S(CC ) .GT .0.1D-20)     ISIG 

10  S( IK)     =    CC 

11  CONTINUE 

IF( ISIG. EC. 1)     GO    TO    15 


=    0 


IF    THE     INVERSE     IS    STILL    NOT    GGOD    ENOUGH    THE    NEW    TRIAL     IN- 
VERSE    IS    PJT     IN    H    AMD    THE     PROCESS     IS    REPEATED 

DO    13    K=2,N 

KK    =    K    -    1 

DO    12    I=1,KK 

I  K  •  =    K    +    N*  (  I     -    1 ) 

12  H(  IK  J     =     S(  IK) 

13  CONTINUE 
GO  TO  14 

15  W  R  I  T  E  (  6 ,  1 0  0 1 )  KIT 
RETURN 

16  WRITE(6, 1000) 
STOP 

1000  FORMAKlOX.'SrOP.  I ER  =- 1 ,  AND  S  CANNOT  BE  FACTORED.') 

1001  FORMAT  I •-• tlOXt1 INVERSE  FOUND  AFTE R»  , I  2 , • I TE RATI ONS . «  ) 
END 


SUBROUT INE  OS  I NV ( A  ,  C  » N  ,  EPS , I ER ) 


DSINV  INVERTS  A  GIVE\  SYMMET 
(MCCIFIED  IN  THE  PRESENT  C3N 
OF  THE  UPPE3  TRIANGJLAR  FACT 
DESCRIPTION!  OF  PARA 
A:  D0U3LE  PRECI 
GIVEN  SYM 
COEFFICIE 
THE  RESULTAN 
D0U3LE  PRECI 
C:  MATRIX  IN  WH 
OF  A  IS  STOR 
PROGRAM 
N:  THE  NUMBER  0 

MATR  IX 
EPS:  SIMGLE  PRE 
USED  AS  A 
LOSS  OF  SI 
IER:  RESULTING 
IER=0 
IER=-1 


RIC  POSITIVE  DEFINITE  MA 

TEXT  TO  PROVIDE  ONLY  THE 

OR  OF  A) 

METERS : 

SIGN  UPPER  TRIANGULAR  PA 

RIC  POSITIVE  DEFINITE  N 

MATRIX.   ON  RETURN   A  CO 

T  UPPER  TRIANGULAR  MATRI 

SION 

ICH  THE  UPPER  TRIANGULAR 

ED  FOR  LATER  USE  IN  THE 


TRIX 
INVERSE 


RT  OF 
BY  N 
NT A  INS 
X  IN 

FACTOR 
CALLING 


F  ROWS  (COLUMNS)  IN  THE  GIVEN 

CI  SIGN  INPUT  CONSTANT  WHICH  IS 

RELATIVE  TOLERANCE  FOR  TEST  ON 

GNIFICANCE 

ERROR  PARAMETER  CODED  AS  FOLLOWS 

NO  ERROR 
NO  RESULT  BECAUSE  OF  WRONG  INPUT 
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PARAMETER    N    OR    BECAUSE    SOME     RAD- 
ICAND    IS    NON-POSITIVE    (MATRIX    A 
IS    NOT    POSITIVE    DEFINITE,     POS- 
SIBLY   DUE     TC    LOSS    OF     SIGNIFI- 
CANCE) 
IER=K       WARNING    WHICH    INDICATES    LOSS    OF 
SIGNIFICANCE.       THE    RADICAND 
FORMED    AT    FACTORIZATION    STEP    K+l 
WAS    STILL    POSITIVE     BUT    NO    LONGER 
>     |  EPS- A (K+l, K+l)  | 


DOUBLE     PRECISION    A,  DIN,  WORK, C 
DIMENSION    A(l)  ,C(1  ) 

FIND    THE    UPPER    TRIANGULAR    FACTOR    OF    A 
CALL     DMFSD(A,\J  ,EPS,IER) 
IF(IER)     9,1,1 

INVERT    UPPER    TRIANGJLAR    MATRIX    T 
PREPARE     INVERSION    LOOP 

1  IPIV=N*  (N+l )/? 
1ND=IPI V 

SAVE    T    BY    DUPLICATING     IT     IN    C 
DO    7     1  =  1,  I  P  IV 
7    C( I)     =     A( I ) 

ITITIALIZE     INVERSION-LOOP 
DO    6     1  =  1  .N 

DIN=1.D0/A ( I  PI  V) 
A(  IPI  V)  =DIN 
MIN  =  N 
KEND=I-1 
LANF=N-KEND 
IF(KEND)     5,5,? 

2  J=IND 

DO    4    K=1,KEND 

INITIALIZE    ROW-LOOP 
WORK=O.DO 
MIN=MIN-1 
LHOR=I PIV 
LVER=J 

START     INNER    LOOP 
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DO    3    L  =  LANF,MI"N 

LVER=LVER+1 

LHOR=LHOR+L 

3  WCRK=WCRK+A(Li/ER)*A(LHOR) 
A(  J)  =-l-.'ORK*DIN 

4  J=J-MIN 

5  I  PI  V=IPI  V-MIN 

6  IND=IND-1 
9    RETURN 

END 


SU3R0UTINE  DMF SD { A , N , E PS , I ER ) 

DMFSD  FACTIRS  A  GIVEN  SYMMETRIC,  POSITIVE  DEFINITE  MATRIX 

0  ITS  UPPER  AND  LOWER  TRIANGULAR  FACTORS  (WHICH  ARE  THE 
TRANSPOSE  OF  EACH  OTHER).   THE  SOLUTION  IS  DONE  USING  THE 
IE  ROOT  METHOD  3!=  CHOLESKY. 

DIMENSION  Ad) 

DOUBLE  PRECISION  DPIV,DSUM,A 

TEST  ON  WR3NG  INPUT  PARAMETER  N 
IF(N-l)  12,1,1 

1  IER  =  0 

INITIALIZE  DIAGONAL  LOOP 
KPI V=D 
DO  11  K=1,N 
KPI V  =  KPI V+K 
IND=KP IV 
LEND=K-1 

CALCULATE  TOLERANCE 

T0L=A3S( EPS*SMGL(A(KPIV) ) ) 

START  FACTORIZATION  LOOP  OVER  K-TH  ROW 
DO  11  I=K,N 
DSUM=O.DO 
IF(LEND)  2,4,2 
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START  INNER  LOOP 

2  DO  3  L=1,LEND 
LANF=KPIV-L 
LIND=IND-L 

3  DSUM=DSUM+A{LANF)*A(LIND) 

TRANSFORM    ELEMENT    A(IND) 

4  DSUM=A(IND)-DSUM 
IF(I-K)     10,5,10 

TEST  FOR  NEGATIVE  PIVOT  ELEMENT  AND  FOR  LOSS  OF  SIGNIFI 
CANCE 

5  IF( SNGL(DSUMJ-TOL)     6,6,9 

6  IFIDSUM]     12,12,7 

7  IF(IER)     8,8,9 

8  IER=K-1 

COMPUTE     PIVOT     ELEMENT 

9  DPIV=DSQRT(DSJM) 
At  KPIV] =DPI V 
DPIV=1 .DO/DP W 
GO    TO    11 

CALCULATE    TERMS    IN    ROW 

10  A(IND)=DSUM*DPIV 

11  IND=IND+I 

RETURN 

12  IFR=-1 
RETURN 
END 
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SUBROUTINE    T MPRD ( H , U ,C , N t F ) 

TMPRD    FORMS    THE    MATRIX    PRODUCT    U**(-T)*H*U**(-1)     AND    THEN 
MULTIPLIES    IT    BY    THE     SAME    FACTOR    AS    WAS    USED    TO    SCALE    S 

IMPLICIT    REALMS     (A-HT3-Z) 

DIMENSION    H(l)  ,U(1)  ,C(  1) 

LIM    =    N*(N    +    1  )/2 

DO    3    L=1,N 

DO    3    K=L,N 

LK    =    K*CK    -    l)/2    +    L 

UHU    =    O.ODO 

DO    2    I=1,L 

HU    =    O.ODO 

IL    =    L*(L    -    I) 12    +     I 

DO    1    J=1,K 

U    =    J*(J    -    1J /2    +     I 

IFCI.GT.JJ     IJ    =     I* ( I    -    l)/2    +    J 

J  K    =    K  *  ( K    -    1 )  /  2    +    J 

1  HU    =    HU    +    H(IJ)*U(  JKJ 

2  UHU    =    UHU    +    U(  IL)*HU 

3  C(LKJ     =    UHU 
DO    4     1=1, LIM 

4  H( I )    =    C(I)*F 
RETURN 

END 
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SUBROUTINE  DJACVT  (  A ,N , NOYE S , E I VU, E I VR ,NDI M) 


DJACVT  CALCUL 
OF  A  REAL,  SY 
THRESHOLD  MET 
ALSO  CALCULAT 
POOLE,  COMPUT 
DESC 


ATES 
MMETR 
HOD. 
ED,  I 
ER  CE 
RIPTI 
A:  TH 
EI 
DI 
N:  OR 
NOYES 


EIVU: 


EIVR: 


NDIM: 


REMARKS 


IN  DOU 

i:  MAT 

CORRE 

F  DESI 

NT  ER, 

ON  OF 

E  REAL 

GENVAL 

MENS  10 

DER  OF 

:  INDI 

NOYE 

NOYE 

THE  E 

TOR, 

LEAST 

NO  PA 

THE  E 

TURNE 

ROW  D 

THE  C 

SION 

PROGR 

CORRE 

IF  EI 

UNDIM 

THE  P 

ME  NT. 

THE  R 

AND  E 

THE  C 


BLE  PRECISION  A 

RIX  BY  USE  OF  T 

S PON  DING  ORTHCG 

RED.   ORIGINAL L 

U.  C,  BERKELEY. 

PARAMETERS 

*8  SYMMETRIC  (SQUARE) 

UES  AND  EIGENVECTORS  A 

N  MUST  3E  NDIM. 

A,  0<N</=NDIM<161 
CATCR  SUPPLIED  BY  USER 
S=0   EIGENVECTORS  ARE 
S=l   EIGENVECTORS  ARE 
IGENVALUES  ARE  RETURNE 
WHICH  MUST  BE  DIMENSIO 

N  IN  CALLING  PROGRAM. 
RTICULAR  ORDER. 
IGENVECTORS  (IF  REOUES 
D  AS  THE  COLUMNS  OF  TH 
IMENSION  OF  WHICH  MUST 
ALLING  PROGRAM.  THE  C 
MUST  BE  AT  LEAST  N  IN 
AM.  THE  VECTOR  EIVR (I 
SPONDS  TO  THE  EIGENVAL 
GENVECTORS  ARE  NOT  REQ 
ENSICNED  VARIABLE  SHCU 
CSITICN  OF  EIVR  IN  THE 


LL  THE 

HE  JACC 

CNAL  EI 

Y  PROGRAMMED  BY 


EIGENVALU 

BI  VARIAB 
GENVECTCR 


MATRIX  WH 
RE  SOUGHT 


ES 
LE 
S  ARE 

W. 


OSE 
.  ROW 


NOT  WANTE 
WANTED 
D  IN  THIS 
NED  TO  AT 
THEY  AR 

TED)  ARE 
IS  MATRIX 

BE  NDIM 
OLUMN  DIM 
THE  CALLI 
» J) t  1=1. 
UE  EIVU(J 
JESTED,  A 
LD  BE  PUT 

CALL  STA 


CW  DIMENSION  (FIRST  DIMENSION) 
IVR  IN  THE  DIMENSION  STATEMENT 
ALLING  PROGRAM.   NDiM<161 


D 

VEC- 

E  IN 

RE- 

,  THE 

IN 

EN- 

NG 

N 

). 

N 

IN 
TE- 

OF  A 
OF 


1) 

2) 

3) 


4) 


REFERENCES 

VDN 


ATRIX  A  IS  DESTROYED  DURING  CALCULATION. 

LL  EIGENVALUES  ARE  CALCULATED  EACH  TIME  DJACVT 

S  CALLED 

MESSAGES  MAY  BE  SENT  TO  STAN- 
BY  DJACVT.   THESE  INDICATE 
AND  ARE  GENERALLY  SELF- 


ARICUS  DIAGNOSTIC 
ARD  OUTPUT  MED  I 
MPROPE*  ARGUMENTS 
XPLANATORY. 
IMI7  Cr  160  FOf 
HANGING  STATEMENT 


(AND 
21 


NDMI)  COULD  BE  RAISED  BY 


NEUMAN,  J.,  AND  GOLDSTEIN,  "THE  JACOBI  METHOD 
FOR  REAL  SYMMETRIC  MATRICES",  JOURNAL  OF  ACM,1, 
1959 
WILKINSON,  J.  H.,  "THE  ALGEBRAIC  EIGENVALUE  PROBLEM 
OXFORD  UNIVERSITY  PRESS,  1965,  P.  277  FF. 

IMPLICIT  REAL* 8 ( A-H,0-Z) 

DIMENSION  A( NDIM, NDIM) , E I VU ( NDI M) , E I VR ( N DI M, NDIM ) 

IF(N-1 ) 20, 23, 2  1 

20  PRINT  22, N 

22  FORMAT («   N=  '.IS,1  IS  TOO  SMALL.  LIMIT  IS  1.  RETURN 
1T0  CALLING  ROUTINE  FROM  DJACVT. •) 

RETURN 

23  PRINT  24,  A( 1, 1) 

24  FORMAT  ( 24H  IN  DJACVT  ,  MATRIX  A  =  ,E14.6) 
RETURN 
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21  IF(N-160)1,1,3 

3  PRINT  2»N 

2  FORMAT! •   N=  ',13,'  IS  TOO  LARGE.  LIMIT  IS  160.  RETURN 
1T0  CALLING  ROJTINE  FROM  DJACVT.') 

RETURN 

1     IF{NOYES)99,102 ,99 

99    CONTINUE 

DO    101     J=1,N 

DO    100     1=1, N 

100  EIVR( I , J )=0.0 

101  EI VR( J, J)=1.0 

102  ATOP=0. 
D0'll2    J  =  1,N 
DO    111     1=1 , J 

IF(A(  I, J  )-A(J,  I  )  )90,103,90 
90    PRINT    106, N,N 
106    FORMAT (14H     IN     DJACVT    ( A (  , I  3 , 1H , , I  3 , 3H ) )  , ) 
PRINT     10  8,1 , J, J, I 

108  FORMAT!'     A ( «  , I  3  ,  ' ,  ' ,  13 ,  '  )     AND    A (  • ,  I  3,  • ,  •  ,  I  3,  '  )     WERE 
1UNEQUAL.     TNEY    WERE     REPLACED    BY    THEIR    MEAN     IN    DJACVT') 

A(  I  ,J)=.5*(A< I  , J )+A(J,  I )  ) 

A(  J,  I)=A( I  ,J) 

103  CONTINUE 

IF(ATOP-DABS(A( 1, J) ) ) 104,111,111 

104  AT0P=DA5S( A( I , J) ) 

111  CONTINUE 

112  EIVU( J )=A{ J ,J) 

IF (ATOP) 109,139,113 

109  PRINT  110 

110  FORMAT! 26H  IN  DJACVT,  MATRIX  A   =  0  ) 
RETURN 

113  AVGF  =  DFL0AT(N*(N-1)  )*.  55 
D=0.0 

DO  114  JJ=2,N 

DO  114  11=2, JJ 

S  =  A(  II-l, JJ)/ATOP 

114  D=S*S+D 
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DSTOP=(1.D-06) *D 

THRSH    =DSQRT    (  0/ AV GF) *  ATOP 

115  IFLAG=0 

DO    130    JCCL=2,N 

JC0L1=JC0L-1 

DO    130    IR0W=1,JC0L1 

AIJ  =  A(  IROWtJCOU 

IF(DABS(AI J) -THRSH) 130,130,117 

117  AI I=A( I  ROW, IROW) 
AJJ=A( JCOL,JCDL) 
S=AJJ-AI I 

IF(DABS( AIJ )-l .  D-09*-DABS(S) ) 130, 130, 118 

118  IFLAG=1 

IF{ l.D-10*DABS (AIJ)-DABS(S))116,119,119 

119  S=. 7071067811865475 
C  =  S 

GO    TO    120 

116  T=AIJ/S 
S=0.25/DSQRT(0.25+T*T) 
C=DSQRT(0. 5+S) 
S=2.*T*S/C 

120  DO    121    I=1,IR3W 
T  =  A(I  , I  ROW) 
U=A( I, J COL ) 

A{ I ,IROW)=C*T-S*U 

121  ACIfJCOL)=S*T+C*U 
I2=IR0W+2 

IF (I2-JCCL) 127,127,123 
127    CONTINUE 

DO    122    1=12, J:OL 

T  =  A(  1-1 , JCOL) 

U  =  A(  IRCW,I-1) 

A(  1-1, JCOL)=S*U+C*T 

122  ACIR0W.I-1)=C*U-S*T 

123  AC  JCOLtJCOU=S*AIJ+C*AJJ 

A( IROW, IROW)=C*A( IROW,  IROW )-S* ( C*AIJ-S*AJ J ) 
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DO    124    J=JCOL,N 
T=A(IROW,J) 

U=A( JCCL,J ) 
ACIROWtJ)=C*T-S*U 

124  AUCCL,  J)  =  S*T+C*U 
IF(NOYES) 131,126,131 

131  CONTINUE 

DO  125  1=1, N 

T=EIVR(  I  ,IR01  ) 

EIVR( I , IROW)=C*T-EIVR< I  ,JCOL)*S 

125  EIVRCI  ,  JCCL)=S*T  +  E I VR(I  ,JCOL)*C 

126  CONTINUE 
S=AI J/ATOP 
D=D-S*S 

IF( D-D  STOP) 126  0,12  9, 12  9 
1263  D=0. 

DO  128  JJ=2,N 
DO  128  11=2, JJ 
S=A( I I-1,JJ)/AT0P 

128  D=S*S+D 
DSTOP=(1.D-06) *D 

129  THRSH=DSCRT(D/AVGF)*ATOP 

130  CONTINUE 

IFf IFLAG1115, 134,115 
134  T=A( 1,1) 

A(1,1)=EIVU(1) 

EIVUC1)=T 

DO  132  J=2,N 

T=A( J, J ) 

A( J, J)=EIVU( J) 

EI VU( J) =T 

DO  132  1=2, J 

132  A(I-1, J)=A( J,I-1) 

133  RETURN 
END 
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SUBROUTINE  MPRD ( A, B , C, N ) 


MPRD  FORMS  THE  MATRIX  PRODUCT  A*B.   THE  RESULT  IS  RETURNED 
IN  B.   C  IS  USED  FOR  TEMPORARY  STORAGE. 

IMPLICIT  REALMS  (A-H,0-Z) 

DIMENSION  At  1)  ,6(1)  ,C(  1) 

DO  2  I =1 ,N 


1 


DO  2  K=1,N 

IK  =  N*(K-1) 

CC  =  O.ODO 

DO  1  J=1,N 

IJ.  =  N*  ( J-l)  +  I 

JK  =  N*(K-1)  *  J 

1  CC  =  CC  +  At IJ)*B(JK) 

2  CCIK1  =  CC 
LIM  =  N*N 

DO  3  1=1, LIM 

3  B(  I)  =  Ct  I  ) 
RETURN 

END 
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SUBROUTINE  INF  256 ( FCT , F ) 

INT256  IS  A  256-POIVJT  GAUSS-LEGENDRE  QUADRATURE  USING  AN  EX- 
TERNALLY DEFINED  FUNCTION.   THE  SAMPLE  POINTS,  WEIGHT  VALUES 
AND  INTEGRATION  SEGMENTS  ARE  READ  INTO  THE  CALLING  PROGRAM 
FROM  AN  EXTERNAL  SOJRCE 

IMPLICIT  REAL*8  (A-H.O-Z) 

CCMM0N/L256/X( I2S) , A(128 ) 

COMMON/ CD/ DP C( 15), DMC( 15) 

F  =  0.  ODO 

DO  100  1=1 ,128 

ZZZ  =  O.ODO 

XX  =  X( I ) 

AA  =  A( I ) 

TO  MINIMIZE  COMPUTATION  TIME  THE  INTEGRAND  IS  COMPUTED  AT 
THE  SAME  RELATIVE  POINT  IN  ALL  15  INTEGRATION  INTERVALS 
FIRST  AND  THEN  THE  SUM  OF  THESE  VALUES  IS  MULTIPLIED  6Y  THE 
WEIGHT  FACTOR  FOR  THAT  POINT 

DO  10  J =1,15 

THE  DPC  AND  DMC  ARE,  RESPECTIVELY,  (D+CJ/2  AND  (D-C)/2, 
WHERE  D  AMD  C  ARE  THE  UPPER  AND  LOWER  INTEGRATION  LIMITS 

DP  =  DPC(J) 

DM  =  DMC (J  ) 

CHANGE    THE    VARIABLE    TO    CONFORM    TO    THE    -1,+1     INTERVAL    RE- 
QUIRED   IN    GAUSS-LEGENDRE    QUADRATURE 

P    =    XX*DM 

ZP  =  DP  +  P 
ZM  =  DP  -  P 
10  ZZZ  =  ZZZ  +  DM* (FCT  (ZP  )  +  FCT(ZM)) 
100  F  =  F  +  ZZZ*AA 

RETURN 
END 


DOUBLE    PRECISION    FUNCTION    CGARG(Z) 

CGARG    EVALUATES    THE     INTEGRAND    AT    THE     DESIGNATED    POINT    FOR 
THE    GAMMA    FUNCTION    INTEGRATIONS 

IMPLICIT    REAL*8     (A-H,0-ZJ 
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COMMON/CONST  /  P'H  1,PH  2,  FRl,FR2»SCL,TwOP  I  ,FCT2»FCT1  ,A 

ARG    =    A*DLOG(Z) 

CGARG    =     (Z»*3)*DEXP (-Z ) *( FGT 2*DS  IN < ARG )     +    FCTl*DCOS(AR 

RETURN 

END 
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SUBROUT INE  TBL2( IKAX.MAXI  , AS, AC, BS,BC) 

TBL2  TABULATES  THE  REQUIRED  BOUND-FREE  INTEGRALS 
IMPLICIT  REAL*8  (A-H,0-Z) 

COMMON/ CCN ST/ P HI, PH2,FR1,FR2, ASCL , P I , C A , CT , AAA 
COMMON/WAVP/WAV 
COMMON/TRM/ALFA, ETA 
COMMON/ ATNC/Z 
COMMON/RHO/N 

DIMENSION  ACQ  }  ,  AS  (1),BC(1),  BS(  1) 
DIMENSION  Al( 144) , A2( 144) 
DIMENSION  SI  (13)  , C2 ( 13 ) , SS ( 1 1 ) , CC ( 1 1  ) 
DIMENSION  X01500) 
ALFA2  =  ALFA  +  Z 

CALCULATE  THE  FREOJENCY  AND  PHASE  FACTORS  IN  THE  SINUSOIDAL 
TERMS 

FR1  =  WAV/ALFA2 

PHI  =  AAA*DLOG( .2D1*FR1 )  +  ETA 

FR2  =  . 2D1-WAV/ALFA 

PH2  =  AAA-'DLOGl  .2D1*FR2)  +  ETA 

ASCL  =  0.5DO*ALFA/ALFA2 

WRITE (6, 3000)  PHI , FR1 , PH2 , FR2 , ASCL 

IF  Z=l  THE  NUMERICAL  INTEGRATIONS  NEED  NOT  BE  DONE 
IF( AAA. EO. 0.000)  GO  TO  104 

FIND  THE  ZEROES  OF  THE  INTEGRAND  FOR  II 
CALL  ZER0(1,L, XO, 10) 

EVALUATE  THE  II  INTEGRALS 
DO  100  I=1,IM*X 
N  =  I  -  2 

CALL  CS INT ( 1,L, IQ,XO,ZZ, ASCL ) 
1  JO  SKI)  =  ZZ 

FIND  THE  ZEROES  OF  THE  INTEGRAND  FOR  IS 
CALL  ZERG(2,L, XO, IQ) 
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EVALUATE  THE  IS  INTEGRALS 
DC  131  1=1, MAX  I 
N  =  I  -  2 
CALL  CSINT(2,L,I 0 , XO, ZZ , 0.  1 Dl ) 

ioi  ssm  =  zz 

FIND  THE  ZEROES  OF  THE  INTEGRAND  FOR  13 
CALL  ZERO( 3,L, XO, IQ) 

EVALUATE  THE  13  INTEGRALS 
DO  102  1=1 , IMAX 
N  =  I  -  2 
CALL  CSINT (3,L, IQ,XO, ZZ,ASCL ) 

102  C2( I )  =  ZZ 

FIND  THE  ZEROES  Oc  THE  INTEGRAND  FOR  IC 
CALL  ZER0(4,L,XJ),  IQ) 

EVALUATE  THE  IC  INTEGRALS 
DO  103  I=1,MAXI 
N  =  I  -  2 
CALL  CSINT{4,L,I0,X0,ZZ,0.1D1J 

103  CC( I )  =  ZZ 
GO  TO  105 

FOR    Z=l     EVALUATE    THE    REQUIRED    INTEGRALS    EXACTLY 

104  CALL    HINTl S1,C2, SS,CC, IMAX,MAXI ) 

105  WRITE(6,2001 ) 

WRITE(6,20  00)(J,S1(J),C2(J),J=1,IMAX) 
WRITE(6, 2002) 

WRITE{6,2000)(  J,SS  (J), CC(J ) , J=1,MAXI ) 
MAXII     =    MAXI     +     1 

INITIALIZE  THE  AS,  AC,  BS  AND  BC  ARRAYS 
DO  199  J=l ,1452 
AS(JJ  =  .ODO 
BS(J)  =  .ODD 
AC(JJ  =  .ODO 
199  BC(J)  =  .ODO 
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INITIALIZE  THE  Al  AND  A2  ARRAYS 
DO  200  J=l,144 

A1UJ  =  .ODO 

200  A2(J)  =  .ODO 

IMAX    =    MAXI 
LLMAX    =     MAXI 
NNM    =    NNMAX    +     1 
LLM    =    2* LLMAX    -    1 
AP    =    0.5D0-AL-A 
ALFAP     =    0. 101/ AP 
ALFA4    =    0.2D1/ (0.2D1*Z    +    ALFA) 

USING    THE    IS    AND    IC     INTEGRALS    FOUND    ABOVE    EVALUATE    AS(NU, 
LAMDA,2)     AND    AC( NU , I AMDA , 2 ) 

DO    201    NN=1, NNMAX 

SSN    =    SS(NN) 

CCN    =    CC(NN) 

AP    =    AP*ALFAP 

F    =    0.1D1 

ALFA4P     =    0.1D1 

DO    2CL     LL=1»LLMAX 

ALFA4P    =    ALFA4P*ALFA4 

F    =    F*DFLOAT{LL) 

IFHNN+LL.LT.'* ) .OR . (NN+LL.GT . IMAX ) )    GO    TG    201 

INA    =    NN    +    NNM*(LL    +    LLM) 

TAA    =     F^ALFA4P*AP 

AS(INA)       =TAA*SSN 

AC(INA)        =TAA*CCN 

201  CONTINUE 
N1MAX    =    MAXI I 
LI  MAX    =    MAXI I 
IMAX I    =     IMAX    f     1 

USING    THE     II    AND    13     INTEGRALS    FOUND    ABOVE    EVALUATE    AKNU, 
LAM  DA* 2)     AND    A2 ( NU , L AMDA , 2 ) 

DO    202    NN=1,N1MAX 

DO    232     LL=1 ,L1 MAX 

IF( (NN+LL.GT. IMAXI ) .OR .(NN+LL.LT. 5) )     GO    TO    202 
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IN  12    =    NN    +    N1MAX* ( LL    -    1) 
CALL    A12SUM(NN,LL,AA,S1) 

Alt  INI12)        =    AA 

CALL    A12SUM(NN,LL,AA,C2) 

A2( IN12)        =    AA 

202  CONTINUE 
NNMAX    =    MAX  I  I 
LLM1    =    LLMAX    -     1 

NLM    =    NNMAX* LLMAX    -    1 

NNM    =    NNMAX    -    1 

LNM    =    NNMAX* (_LMAX    +    1)     -    2 

USING    ThE    RECURSION    RELATIONS    EVALUATE    AS ( NU , LAMDA, 1 ) i 
AC(NLU  LAMDA, 1),     E>S  (  M  U  ,  LAMDA  ,  1 )     AND    BC  (  NU  t  LAMDA  ,  1  ) 

DO    2  33    NN=2f NNMAX 

DO    2  03     LL=1, LLMAX 

IFC(NN+LL.GT.IMAXI J .OR . (NN+LL. LT.5) J    CO    TO    203 

INI    =    NN    +    NNMAX* (LL    +    LLM1J 

INA1    =     INI    +    NLM 

INLA    =    LL   +    N1MAX* ( NN    -    1) 

.A    =    LL     +    I     +    N1MAX*(NN    -    2) 

AS(I\1)       =       AS(INAl)     +       AKINLA)    -       AHINNA) 

ACIIN1)        =       AC(INAl)     +       A2(INLA)    -       A2UNNA) 

IF((NN.LT.3).DR.(LL.LT.2) .OR .( NN+LL, LT.6 J . CR.tNN.GT 
1NNM).0R. (LL.GT.LLM) )    GO    TO    20 

INB1    =     INI    +    LNM 

INLB    =    LL    --    1    +    N1MAX*NN 

INNB    =    LL    +    2    +    N1MAX*(NN    -    3) 

BS(ni)       =     (     ASCINB1)     +       AKINLBJ    -       Al  (  I  NNB)  ) /.  3  Dl 

BC(INl)       =     (    ACCINB1J    +       A2UNLB)    -       A2(  INNB)  J/.3D1 

203  CONT  INUE 
MMMAX    =    MAXI 
MAXI3    =     IMAX    *-    3 
MMM    =    MMMAX    -    2 

NN1    =    2*(NNMAX*LLMAX    -    1) 
NN2    =    2*NNMAX*LLM1 
NN3    =    NNMAX-LLM    -     1 
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USING    THE     RECURSION    RELATIONS    EVALUATE    AS ( NU , LAMDA, MU ) , 
AC(NU,  LAMDA, MU),     BS ( NU , L AMDA , MU)     AND    BC { NU , LAMDA , MU)     FOR 
MU    >    2 

DO    204    MM=4tMMMAX 

FM31    =    DFLOAKMM    -    3)/DFL0AT(MM    +    1) 

DO    204    NN=2fNNMAX 

DO    204    LL=1,LLMAX 

IF ( (NN+LL+MM.GT.MAXI3) .OR. {NN+LL.LT.4)  )     GO    TO    234 

IN    =    NN    +    NNMAX*((LL    -     1)     +    LLMAX*(MM    -     1)) 

INN    =    IN    -    NN1 

INL    =     IN    -    NN2 

INNL    =     IN    -    HH3 

AS(IN)        =       AS(INN)     +       AS(INL)    -       BS (  INNL ) * . 2D  1 

AC(IN)        =       AC(INN)     +       AC(INL)     -       BC( INNL )* .2D1 

IF{  (  LL.  LT.  2)  .DR.  (NN.LT.3)  .OR. ( LL  +  MM.LT.6)  .OR. (MM.GT. 
1MMM) )     GO    TO    20  4 


BS(IN)        =    (     BS(INN)  +       BS(INL) 

BCIIN)       =     (     BC(INN)  +       BC(INL) 
204    CONTINUE 
RETURN 
2000    F0RMAT(20X,  I2t 2D30.  16) 

2JJ1    FORMATt ////,34X»«  S1',24X,« 

2002     FORMAT (////, 3^X, •  SS',24X,' 


AS ( INNL)*.2D1)*FM31 
AC( INNL )*.2D1)*FM31 


C2' ,// ) 
CC  ,//) 


3000  FORMAT! »-• ,10X,»PH1  =  '  ,  P23. 16 , / , « 0'  , 1  OX , ' FR1 
116,/ . « 0' ,10X, • PH2  =  • ,D23. 16,/ , »0» , lOXt ' FR2  = 
2,/,  '  J»  ,1 JX,' ASCL  =  ',022.16,/,'-') 

END 


=  ' ,D23. 
• ,D23.16 
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SUBROUTINE    ZER  0(  I  C  S  ,  J  ,  XO  ,  I  Q) 

ZERO    FINDS    ALL    THE    ZEROES    OF    SIN(A*X     +    B*LN(X)     +    C)     OR 
COS(A*X    +    B*LN(X)     +    C)     BETWEEN    X  =  XSM    AND    X=SLM 

IMPLICIT    REAL*8     (A-H,0-Z) 

COMMON/CCNST/PHl,PH2,FRl,FR2,ASCL,PI ,CA,CI ,B 

COMMON/INT/XSM,FF,BO,XLM 

DIMENS ION    XO(l  ) 

X    =    XSM 

GO    TO    (21,22,21,22 ) , ICS 

21  A    =    FR1 
CC    =    PHI 
GO    TO    1 

22  A    =    FR2 
CC    =    PH2 

1    ARG    =    A*X    +    B<-'DLOG(X)     +    CC 

FIND    THE    SMALLEST     (ALGEBRAICALLY)     INTEGER    MULTIPLE    OF    PI 
SUCH    THAT    A*X    +    B*LM(X]     +    C    >    N*PI 

ARGP     s    ARG/PI 

ITHTA    =     IDINT(ARGP)    -     1 

IT    =     IABS(  ITHTA) 

DETERMINE  WHETHER  THIS  ZERO  STARTS  A  POSITIVE  OR  A  NEGATIVE 
SWING  OF  THE  FUNCTION  AND  RECORD  THIS  FACT 

10  =  MOD C IT, 21 

IF( IQ.EO.O )  13  =  -1 

THTA  =  PI*DFLDAT( I THTA) 

IF  (ICS.GT.2)  THTA  =  THTA  -  0.5D0*PI 

FIND   THE    VALUE    OF    X    CORRESPONDING    TO    THIS     INITIAL    VALUE    OF 
THE    ARGUMEM 

CALL    XZERO(THT A,X, A,B,CC) 

X0(1)     =    X 

REPEAT  THIS  PROCESS  UNTIL  X  >  XLM  OR  UNTIL  500  ZEROES  HAVE 
BEEN  FOUND 

DO  2  J=2,500 

THTA  =  THTA  +  PI 

X  =  X  +  PI  /( A  +  B/X) 
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CALL    XZEROCTHT  A,X,  A,B,CC) 

XO(JJ     =    X 

IF    (X.GT.XLM)  GO    TO    3 

2  CONTINUE 
J  =  500 

3  CONTINUE 
RETURN 
END 


SUBROUTINE  XZE RO (THTA , X , A , B ,CC ) 

XZERO  SOLVES  THE  EQUATION  THETA  =  A*X  +  B*LN(X)  +  C  FOR  X  BY 
AN  ITERATION  PROCESS 

IMPLICIT  REALMS  (A-H,0-Z) 

1  XI  =  X*(THTA  -  (CC  +  B*(DLOG(X)  -  1)))/(X*A  +  B) 

IF  THE  LINEAR  APPROXIMATION  RESULTS  IN  A  NEGATIVE  ESTIMATE 
FOR  X  THE  METHOD  FAILS  SO  REPLACE  THE  ESTIMATE  OF  X  BY  ONE- 
.TH  OF  THE  INITIAL  VALUE  OF  X  AND  REPEAT  THE  APPROXIMATION 
DING  A  ■       riMA"      REPEAT  THIS  PROCESS  UNTIL  THE 
ESTIMATE      MES  POSITIVE 

IF(X1 .LE.O.ODD )  XI  =  0.1D0*X 

IF  (DABS(X-Xl)  .LT.O.  1D-12)  GO  TO  2 

X  =  XI 

GO  TG  1 

2  X  =  XI 
RETURN 
END 
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SUBROUTINE  CSI\T(ICSfM,IQ,XO,F,P) 

CSINT  CARRIES  OUT  THE  NUMERICAL  INTEGRATION  OF  THE  BOUND- 
FREE  INTEGRALS  FOR  I    >  1 

IMPLICIT  REAL*8  (A-H,0-Z) 

EXTERNAL  FQN, FCN , DS IN, DCOS 

COMMON/KONE/C, D,FFA 

COMMON/ I NT/XSM,FF, B0,XLM 

COMMON/ GEZRO/F  A,  B,  AMP,  ICCIQI 

COMMON/RHO/N 

DIMENSION  STS( ^)  ,TI  (3)  , AC(4) , AD (4) , AES(4) , XO(M) 

ICC  =  ICS 

TS  =  O.ODO 

SINT  =  O.ODO 

AES(l)  =  Q.1D1 

DO  4  LL=1,3 

4  AES(LL+1)  =  AESILL)  +  P 
FN  =  DFLOAT(N) 

SINCE  THE  REGION  WHERE  THE  INTEGRAND  IS  SIGNIFICANTLY  DIFF- 

ENT  FROM  ZERO  IS  DIFFERENT  FOR  DIFFERENT  N,  ADJUST  THE 
PC'i  iERE  THE  CAREFUL  NUMERICAL  INTEGRATION  STARTS  AND 
ENDS  TO  TAKE  ADVANTAGE  OF  THIS 

XS  =  XSM 

IF(N.GT.l)  XS  =      FF*(FN  -  D . 1 Dl ) *( FN** J  .6D 3 ) 

XL  =  0.25D1*FM  +  BO 

IF(N.LE.  1)  GO  TO  7 

DO  5  J  =  l ,M 

IF(XO( J )  .GE.XS)  GO  TO  6 

5  CONTINUE 

6  JJ  =  J  -  1 

IF( JJ.EQ.O)  JJ  =  1 

IOI  DETERMINES  WHETHER  THE  FUNCTION  IS  STARTING  A  POSITIVE 
OR  NEGATIVE  SWING  OM  THE  NEXT  HALF-CYCLE 

IOI  =  MOD( JJ,2 ) 

IF( IQI .EO.l)  IOI  =  IQ 

IF( IOI.EO. 0)  I  01  =  -IQ 

GO  TO  8 
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7  JJ    =    1 
101    =    10 

8  D    =    X0( JJ) 

IF    (N.EQ.-l)     30    TO    99 
DO    9    LL=1,4 

IF(D.EO.O.ODO)     GO    TO    100 
AD(LL)     =    DEXP{-AES  (LL)*D)/D 
GO    TO    9 
130    AD(LL)     =    O.ODO 

9  CONTINUE 

99    DO    10    J=1,M 

J I     =    M    +    1     -    J 

IF    [XO(JI).LT.XL)    GO    TO    11 

10  CONTINUE 

11  IF(JJ.GE.JI)    JI    =    M    -     1 
DO    17    K  =  JJ  ,JI 

C    AND    D    REPRESENT    SJCCESSIVE     ZEROES    OF    THE    TRUE     INTEGRAND 
C    =    D 
D    =    X0(K  +  1  ) 

FA     IS    THE    FREQUENCY    OF    THE    SIN    FUNCTION    WITH    ZEROES    AT    C    AND 
D 

FA    =    0.314159Z6D3589793D1/ (D    -    C) 

FFA    =     FA 

DIO  =  DFLCATU  01  ) 

IF  THE  TRUE  INTEGRAND  IS  POSITIVE  THE  APPROXIMATE  FUNCTION 
IS  REDUCED  IN  AMPLITUDE  BY  10?  AND  IF  NEGATIVE,  INCREASED 
BY  10% 

AMP  =  0.  1D1  +  0. 1D0  ^DIC 

IF(N.EO.-l)  GD  TO  14 

DO  12  LL=1»4 

AC( LL)  =  AD(LL) 

12  AD(LL)  =  DEXP(-AES(LL)*D)/D 

FIRST  THE  INTEGRAL  OF  THE  APPROXIMATE  FUNCTION  IS  EVALUATED 
DO  13  LL=1,4 
AE  =  -AES( LL) 
AEC  =  AC(LL  ) 
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AED  =  AD(LL) 

CALL  CSSUMCDtC ,FA, TERM, AE,AEC, AED) 
13  STS(LL)  =  TERM 


TS  =  TS  +      DI0*AMP*(STS(1) 
1  -  STS(3)J  -  STS(4)  J 

GO  TO  16 

14  UL  =  0.1D1 


0.3D1* (STS(2) 


FOR  N=-l  THE  APPROXIMATE  FUNCTION  MUST  ALSO  BE  EVALUATED 
NUMERICALLY 

DO  15  LL=1,3 

BL  =  UL 

UL  =  BL  +  P 

CALL    D0G32(BLtUL tFQNtTT) 

15    TI(LL)     =    TT 

TS    =    TS    -  DIQ*AMP*FA*(TI ( 1)     -    0.2D1*T1(2)     +    TI13I) 


NEXT  EVALUATE  THE  INTEGRAL  OF  THE  DIFFERENCE  BETWEEN  THE 
TRUE  INTEGRAND  AND  THE  APPROXIMATE  ONE  NUMERICALLY  USING  A 
32-POINT  GAUSS-LEGENDRE  QUADRATURE 


16  DD  =  DFLOATd  +  (1 
CC  =  DD  +  0.101 


IOD/2  ) 


B    IS    THE    PHASE    OF    THE     APPROXIMATE     FUNCTION,     SO    CHOSEN    TO 
MAKE    THE    FJNCTION    BEHAVE     PROPERLY 

B    =    FA*(DD*D    -    CC^C) 

CALL    D0G32(  CD,  FCN,  FG) 

SINT    =    SINT    +    FG 

17    IOI    =    -IQI 

FINALLY    THE    REGION    PROM    0    TO    XS    AND    FROM    XL    TO     100     IS     INT- 
FGPATED    WITHOUT    USIM3     THE    APPROXIMATE     FUNCTION,     UTILIZING    A 
512-PC1NT    GAUSS-LEGENDRE    QUADRATURE 

IF    (  ICS.GT.2)    GO    TO     18 

D    =    X0( JJ) 

IF(D.EQ.O.ODO)     GO    TO    24 

CALL    DC512( O.ODO,D ,DSIN,ZI ) 

GO    TO    25 

24  ZI    =    0.00  0 

25  C    =    XOCJI+IJ 
IF(C.GE.0.6D2)     GO    TO    20 
CALL    DQ512(C,0.6D2 ,DSIN,UI ) 
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CL    =    0.6D2 
GO    TO    21 

20  CL    =    C 

UI     =    O.ODO 
IF(C.GE.0.1D3)     GO    TO    26 

21  CALL    DQ512(CL,0.1D3,DSIN,UJ) 
GO    TO    19 

18  D    =    X0( JJ) 

IFtD.EO. O.ODO)     GO    TO    28 
CALL    DQ512(0.DD0tD, DCOS, ZI ) 
GO    TO    29 

28  ZI     =    O.ODO 

29  C    =    XOUI  +  1) 
IF(C.GE.0.6D2)     GO    TO    22 
CALL    D0512(Ct3 .6D2, DCOS,UI ) 
CL    =    0.6D2 

GO    TO    23 

22  CL    =    C 

UI     =    O.ODO 
IF(C.GC.0.1D3)     GO    TO    26 

23  CALL    D0512(CL,0.1D3,DC0S,UJ) 
GO    TO    19 

26    UJ    =    O.ODO 

THE     FINAL     INTEGRAL     IS    THE     SUM    OF    ALL     THE     INDIVIDUAL    INTEG- 
RATIONS 

19  F    =    SINT    +    TS    +    ZI     +    UI     +    UJ 
RETURN 

END 
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SUBROUTINE  C SS UM ID , C , A .TERM, AE , AEC , AED) 

CSSUM  FINDS  THE  INTEGRAL  OF  (X**N )* ( EXP (-ALPHA*X ) )*S  IN( AX  +  B ) 
BETWEEN  SUICESSIVE  ZEROES  OF  THE  INTEGRAND 

IMPLICIT  REALMS  (A-H.O-Z) 

COMMON /RHO/N 

LIM    =    N    +     1 

DEN    =    -DSQRT(A*A    +    AE*AE) 

PH    =    DATAN2(A,AE) 

S    =     3.0D0 

F    =    0. 1D1 

DI     =    0. 1D1 

DR    =       0.1 01/ DEN 

PHP    =    PH*DFLC^T(N+2) 

FN    =    D.1D1/DFL0AT(LIM] 

AC    =    AEC 

AD     =    AED 

DO    10    1=1. LIM 

FI     =    DFLCAK  I  } 

FN    =    FN*FI 

DI    =       DI*DEN 

DR    =       DR*DEN 

F    =    F*F 1 

PHP    =    PHP    -    PH 

AC    =    AC*C 

IF ( {AC. EO. J. JD3) .AND. ( I .EQ.l ) )     AC    =    0.1D1 

AD    =    AD*D 

10    S    =    S    +    DSIN(PHP)*FI*DR*(AD    +    AO/F 

TERM    =    FN-S/DI 

RETURN 

END 
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SUBROUTINE    DQ3 32 ( X L , XU , FCT ,  Y) 

D0G32     IS    A    32-PDIM    GA USS-LE GENDRE    QUADRATURE    WHICH    INTEG- 
RATES    AN    EXTERNALLY     DEFINED    FUNCTION    FCT    BETWEEN     THE    LIMITS 
XL    AND    XU 

DESCRIPTION     OF    PARAMETERS 

XL:  LOWER  BOUND  OF  THE  INTERVAL 
XU:  UPPER  BOUND  OF  THE  INTERVAL 
FCT:     NAME    GF    THE    EXTERNALLY    DEFINED    FUNCTION 

TD    BE     INTEGRATED 
Y:    RESULTING    INTEGRAL    VALUE 

DOUBLE    PRECISION    XL ,XU , Y , A, B, C , FCT 

DEFINING  A  AND  B  ENABLES  THE  VARIABLE  TO  BE  CHANGED  T3  CON- 
FORM WITH  THE  -1,  +  1  INTERVAL  REQUIRED  OF  THE  G AUSS-L EGENDRE 
QUADRATURE 

A=..5D0*(  XU+XL) 

B=XU-XL 

C=.49  86  319  3  092  474  078D0*B 

Y=.35093050047350483D-2*  (FCT  (A+O  +  FCT  (A-C ) ) 

C=.49  280  57  5  57  726  3417D0*B 

Y=Y+.  813719736  5452  8350-2*  (FCT  (A+O+FCT  lA-O) 

C=.48  2  38  112779  37  53  2  2D0*B 

Y-Y+.l 26  9  60326 5463 1 C3GD- 1* (FCT (A+C) + FCT (A-C) ) 

C=. 46745303796 886984D0*B 

Y=Y  +  . 17 13693145651 0717D-1*( FCT (A+O+FCT  (A-C J ) 

C=.44816  357788  3  326  36D3*B 

Y=Y+. 21 417949011  1 1 3340D-1* ( FCT [A+O+FCT  (A-C) ) 

C=. 42468380 686 62 8499D0*B 

Y=Y+. 2  5499  32  9S31 188  088  0-1*  (  FCT  (  A  +  O+FCT  (  A-C)  ) 

C=  ,39724189798397120D0*B 

Y  =  Y+. 29  342  3467  3926  77  74  0-1* ( FCT ( A+C ) +FCT ( A-C ) ) 
C=.  365 09 10 5937 0  1448 400 *B 

Y=Y+.  3291111 13  88180923D-1M  FCT  (A+O+FCT  (A-C)  ) 

C=. 3315221334651 376OD0*B 

Y=Y+. 36 1728970 54424253D-1*(  FCT(  A+O+FCT (  A-C)  ) 

C=. 293 8578 7862 0381 16D0*B 

Y=Y  +  .39096947893535153D-1*(  FCT  (A+O+FCT  (A-C)  ) 

C=.2  534499  5446  6114  7  0D0*B 

Y  =  Y+.  4  16559621  1347  3378  0-1*  (  FCT  (A+O+FCT  (A-C)  ) 
C=. 21067563806 53 1767D0*B 
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Y  =  Y+.43  826  J46  5,)2  20190  6D-1* ( F CT ( A+C ) +FCT ( A-C ) ) 

C=.16 59 343011^ 1063  8  2D0*B 

Y= Y+. 45 5 869393 47  881 942 D- I'M  FCT ( A  +  C )+FCT( A-C)  ) 

C=.119  64  36  811ZSJ68  54DJ*B 

Y=Y+.46922199540402283D-1*(FCT( A+C)+FCT( A-C) ) 

C=.722  3  598  079139  82  5D-1*B 

Y=Y+.4  7C19360D3963  743  0D-1*(FCT(A+C)+FCT( A-C) ) 

C=.2415383  2  843  86  915  8D-1*B 

Y=B*(Y+. 4827004425 7363900D-1*< FCT ( A+C ) +FCT ( A-C ) ) ) 

RETURN 

END 


DOUBLE    PRECISION    FUNCTION    FQN(Z) 

FON    EVALUATES    THE    APPROXIMATE    FUNCTION    AT    GIVEN    Z    FOR    RH0=-1 
IMPLICIT    REAL*8     (A-H,0-Z) 
COMMON/MONE/C,  D, FA 

FON    =     (DEXP(-:-Z)     +    DEXP(-D*Z) }/(Z*Z    +    FA-FA) 
RETURN 
END 


DOUBLE    PRECISION    FUNCTION    FCN(Z) 

FCN    EVALUATES    THE    DIFFERENCE     BETWEEN    THE    TRUE     INTEGRAND    AND 
THE    APPROXIMATE    ONE    FOR    GIVEN    Z 

IMPLICIT    REAL*8     (A-H.O-Z) 

COMMON/ CON ST/ P HI,  PH2,FR1,FR2, ASCL,PI ,CA,CI  ,AA 

C0M^1DN/GEZR0/rA,B,F  ,  ICS,ICI 

COMMON/ RHO/N 

GO    TO    (  1,2,  1,2)  ,  ICS 

1  FR  =  FR1 
PH  =  PHI 
P  =  ASCL 
GO    TO    3 

2  FR    =    FR2 
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PH    =    PH2 
P    =    0.1D1 

3  ARGF    =     FR*Z    +    AA*DLOG(Z)     +     PH 
ARGH    =     FA*Z    +    8 

S    =    0.  1D1    -    DEXP(-P*Z) 
ZN    =    0.1D1/  (Z*Z) 
NN    =    N    +    2 
DO    6    I =1 ,NN 
6    ZN    =    ZN*Z 

ENV    =    S*S*S*DEXP(-Z)*ZN 

IF     (  ICS.GT.2)    GO    TO    4 

FCN  =  ENVMDSI  N(  ARGF)  -  F*DS  IN ( ARGH)  ) 

GO  TO  5 

4  FCN  =  ENV* (DCOS( ARGF)  -  F*DSIN( ARGH) ) 

5  RETURN 
END 
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SUBROUTINE    D05 12 ( C , D , FCT , Y ) 

D0512     IS    A    512-POINT    G AUSS-L EGENDRE    QUADRATURE    WHICH    EVAL- 
UATES    THE    BOUND-FREE     INTEGRALS    BETWEEN    THE    LIMITS    C    AND    D. 
THE    EXTERNALLY    DEFIMED    FUNCTION    FCT    IS    EITHER    SIN    OR    COS. 
LOCATION    OF    THE    SAMPLE    POINTS    AND    THE    CORRESPONDING    WEIGHT 
FACTORS    ARE     READ    INTO    THE     MAIN    PROGRAM    FROM    AN    EXTERNAL 
SOURCE  . 

IMPLICIT    REALMS     (A-H,0-Z) 

COMMON/L512/X(256) , A(2  56) 

CCMMON/CO^ST/PHl ,PH2,FR1 ,FR2,ASCL,PI ,CA,CT,AA 

COMMON/GEZRO/FA, B, AMP, ICS, I Q I 

COMMON /RHO/N 

GO.  TO  (  1,2,1,2)  ,  ICS 

1  FR  =  FR1 
PH  =  PHI 
P  =  ASCL 
GO    TO    3 

2  FR  =  FR2 
PH  =  PH2 
P    -    0. 1D1 

3  Y    =     D.JDJ 

THE    VARIABLE    SO    THAT    THE     INTEGRATION    LIMITS    BECOME 
(-1,+1) 

DPC    =    0.5D0-MD    +    C) 

DMC    =     0.5D0*(D    -    C) 

DO    4    I =1,256 

XDC  =   xu)*dm: 

zp   =  dpc  +   xd: 

2M    =    DPC   -    XDC 

ZMN    =    0. 1D1/(ZM*ZM) 

ZPN    =    0.1D1/  UP*ZPJ 
NN    =    N    +    2 
DC    5    J=l ,NN 
ZMN    =    ZMN*ZM 
5    ZPN    =    ZPN*ZP 

SP    =    0.1D1    -    DEXP(-P*ZP) 
SM    =     0.1D1    -    DEXP(-P^ZM) 
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ENVP    =    SP*SP*SP*DEXP(-ZP)*ZPN 
ENVM    =    SM*SM*SM*DEXP(-ZM)*ZMN 
ARGP    =    FR*ZP    *■    AA*DLOG(ZPJ     +    PH 
ARGM    =    FR*ZM    +    AA*DLOG(ZM)    +    PH 
4    Y    =    Y    +    At  I  )*(  ENVP*FCT(ARGP)     +    ENVM*FCT ( ARGM) ) 
Y    =    Y*DMC 
RETURN 
END 
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SUBR3UTINE    HI^J  T(  SI  TC2  ,  SS,CC  , I  MAX t MAX I  ) 

HINT    FVALUATES    THE    BOUND-FREE     INTEGRALS    FOR    THE    CASE    Z=l 

IMPLICIT    REAL*8     (A-H.O-Z) 

EXTERNAL    FACT 

COMMON/ CONST/ PHI , PH2 , FR1 , FR2 , ASCL , P I , CA , C I , AAA 

DIMENSION    RHOl (4), RH02(4) , TriTASK 4) ,THTAC1 (4) , 
1  THTAS2 (4) ,THTSC2(4) 

DIMENSION    SKI  MAX)  ,C2(  IMAX) , SS ( MAXI ) ,CC ( MAXI ) 

Fl    =    FR1**2 

F2    =    FR2~*2 

DO    1    K=l,4 

A    =    DFLOAT (K    -     1 ) 

RHOKK)     =    DSQRTU0.1D1    +    A*ASCL)**2    +    Fl) 

RH02(K)     =    DSQRTU0.1D1    +    A)**2    +    F2 ) 

THTASl(K)     =    DARS  IN ( FR1 / RHO 1 ( K )  ) 

THTACKK)     =    DARCOSC  (  0.  1D1    +    A* ASC L) / RHOl ( K ) ) 

THTAS2(K)     =    DARS  IN(FR2/RH02 (K)  ) 

1    THTAC2(K)    =    DARCOS( ( 0. ID  1    +    A)/RH02(K)J 

SKI)     =    DATAN(FRl)     -    0  .  3 Dl- ( DAT  AN ( FR1/ <0  . 1  Dl     +    ASCL)) 

1  -    DATAN( FR1/(0.1D1    +    0.2D1*ASCL) ) J 

2  -    DATANi FR1/ (0.1D1    +    d.3Dl*ASCL)J 

C2(l)     -    0.15D1*DL0G( ( ( 0.1D1    +    ASCL)**2    +    F1)/((0.1D1 

1  +    0.2D1*  ASCL  )*:-2     +    Fl))     +    J  .  5D0*DL0G  (  (  (D  .  ID  1 

2  +    0.3D1*ASCL)**2    +    FD/C0.1D1    +    Fl)) 

SS(1)     =     DATAN(FR2)     -    0  .  3D1*"  (  DAT  AN  (  FR2/0  .  2D  1  ) 
1  -    DATANC  FR2/0.3DD)    -    DATAN{  FR2/0.4D1 ) 

CC(1)     =    0. 15D1*DL0G(  (0.4D1     +    F2J/(J.9D1    +    F2))    + 
1  0.5D0*DL3G( (0.16D2    +    F2)/t0.1Dl    +     F2 ) ) 

DO    2    J =2, I  MAX 

I     =    J    -     1 

Fl     =    DFLOAT( I ) 

FCT    =    FACT(J    -     2) 

SKJ)     =    FCT*(DSIN(FI*THTAS1  (1)  )  /  RHOl  (1  )**I 

1  -    0.3D1*(DSIN(FI*THTAS1(2) )/RH01(2)**I 

2  -    DSIN(FI*THTAS1(3) )/RH01(3)**I ) 

3  -    DSINt FI*THTAS1(4) ) /RH01(4)**I ) 

C2(J)     =    FCT*(DC0S(FI*THTAC1(1) )/RH01(l)**I 

1  -    0.3D1*(DC0S(  FI~THTAC1(2)  )/RH01(2)**I 

2  -    DC0S(FI*THTAC1(3) )/RH01(3)**I ) 

3  -    DC0S(FI*THTACH4)  )/RH01(4)**I) 

IFU.GT.MAXI)     GO    TO    2 
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SS(J)  =    FCT*(DSIN(FI*THTAS2(1) )/RH02(l)**I 

1  -  0.3Dl*(DSIN(FI*THTAS2(2n/RH02(2)**I 

2  -    DSIN(FI*THTAS2(3) )/RHC2(3)**I ) 

3  -    DSIN(FI*THTAS2(4) )/RH02(4)**I ) 

CC(J)     =     FCT*  (DCOS  (FI-  THTAC2  (1  )  )/RH02(  1)**I 

1  -    0.  3D1-   (DCOS(  FI*THTAC2(2)  J  /RH02(2)**I 

2  -    DC0S(FI*THTAC2(3) )/RH02(3)**I ) 

3  -    DC0S(FI*THTAC214) )  /  RH02  (4  )  •■=- I ) 

2    CONTINUE 
RETURN 
END 
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SUBROUTINE    A  12  SUM  (  N U, L MDA , AA ,VAL) 

A12SUM    EVALUATES    AH  NU,  L  AMDA  ,  2)     OR    A2 ( NU , LAMDA , 2 )     USING    II 
OR     13,     RESPECTIVELY 

IMPLICIT    REAL*8     (A-H,0-Z) 

DIMENSION    VAL(  1  ) 

COMMON /TRM/ ALFA, ETA 

COMMON/ATNO/Z 

ALFA2     =     Z  +    ALFA 

FCTR    =     0.2D1**LFA2/(0.2D1*Z    +    ALFA) 

NU1    =    NU    +     1 
NL    =    NU    +    LMDA 
F    =    0.  1D1/DFL0AT(NU) 
FCT    =    0.1D1/FCTR 
TVAL    =     O.ODO 
DO    10    I =1, NU 
F    =    F*DFL3AT(MU1    -     I) 
FCT    =    FCT*FCTR 
10    TVAL    =    TVAL    +    FCT*VAL(NL-IJ 
AA    =    TVAL/ALFA2**(NL-2) 
RETURN 
END 


•>'  *  i  it  St  ^t  - 
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SUBROUTINE    PHASE  ( AS  ,  AC  ,  BS.  BCMAXI ,  HS,  HCf  E  ) 

PHASE    COMPILES    THE    MATRIX    ELEMENTS    <S I ( H-E )  I  CHI  (  I ) >    AND 
<C  I  (H-E)  I  CHI  (I )> 

IMPLICIT    REAL*8     IA-H,0-Z) 

COMMON/TRM/ALFA.ETA 

COMMON/ATNO/Z 

COMM3N/INDEX/N I . NJ ,LI , LJ ,MI ,MJ ,1 S 

DIMENSION    AS(1  )  ,  AC  ( 1  )  ,  BS  ( 1  )  »  BC  (  1  ) 
DIMENSION    NU(26) ,LMDA{ 26) ,MU( 26) tC0EF(26) 
FN    =    DFLCAT(NJ) 
FL    =    DFLCATtLJ  ) 
FM    =    DFLCATC-'.J  ) 

SET    THE    VALUES    CF     NJ i     LAMDA,    MU    AND    THE    COEFFICIENTS    FOR 
EACH    OF    THE    26    TERMS    OF    THE    MATRIX    ELEMENT 

NU(     1)     =  LJ    +     1 

NU (    3 )     =  NU (  1 ) 

NU(     4)     =  NU(  1) 

NU(     6)     =  NU(1) 

NLU    8J     =  NU(  1) 

NU(     9)      =  NU(  1) 

NU(lv)     =  NU(1) 

LMDA( 12)  =    NU( 1) 

LMDA(26)  =    NU(1) 

NU(     2)     =  LJ 

LMDA( 17)  =    NU( 2) 

NU(21)     =  NU(2) 

NU(     5)     =  LJ    +    2 

NU(20)     =  NU(5) 

NU(22)     =  NU( 5) 

LMDAJ  10)  =    NU(  5) 

LMDAdl  )  =    NU(  5  ) 

LMDAl 14)  =    NU(  5) 

LMDAJ 15)  =    NU( 5) 

LMDA( 16)  =    NU(  5) 

LMDA( 18)  =    NU( 5) 
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LMDA124)  =  NU( 5) 

NU(  7)  =  LJ  -  1 

NUl 10)  =  NJ  +  1 

NUl 12)  =  NUl 10) 

NU( 13)  =  NU( 13) 

NU(15)  =  NU(ID) 

NUl 17)  =  NU( 10) 

NU118)  =  NU(10) 

NU123)  =  NU(  10  ) 

LMDA(  3)  =  NUl 10) 

LMDA122 )  =  NU( 10  ) 

NU.(ll)  =  NJ 

NU{25)  =  NU(ll) 

LMDA(  8)  =  NUl  11  ) 

r;u(  14)  =  nj  +  2 

NU124)  =  NUl  14  ) 

NU(26)  =  NU( 141 

LMDA(  1)  =  NUl 14) 

LMDA(  2)  =  NU( 14) 

LMDAl  5)  =  NUl 14) 

LMDA(  6)  =  NU(14) 

LMDA1  7)  =  NUl 14) 

LMDA1  9)  =  NUl L4) 

LMDA120 )  =  NUl 14) 

NUU6)  =  NJ  -  1 

LMDA1  4)  =  NJ  +  3 

LMDA119J  =  LMDA14) 

LMDAC21J  =  LMDA14) 

LMDA113)  =  LJ  +  3 

LM0AC23J  =  LMDA1  13). 

LMDA125)  =  LMDAU3  ) 

MU1  1)  =  MJ  +  2 

MU1  2)  =  MU(l) 

MU  1  3  )  =  MU  (  1 ) 

MU1  7)  =  MU1 1) 

MU 1  8 )  =  MU 1 1 ) 
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KU{ 10)     =  MU( 1) • 

MU( 11)     =  MU(1) 

MU(12)     =  MU(l) 

MU(16)     =  MU(1) 

MU(17)     =  MU(1) 

MU(     4)     =  MJ 

MU(     5)     =  MU(4) 

MU(    6)     =  MU(4) 

MU( 13)     =  MU(4) 

Ml!  (14)     =  ML)  (4) 

MU(15)     =  MU(4) 

MU119)     =  MIM4) 

MUC20)     =  MU(4) 

MU(21)     =  MU(4) 

MU(22)     =  Ml)  (4) 

MU(23)     =  MUC4) 

MU(24)     =  MU(4) 

MUC25  3    =  MU(4) 

MUC26J     =  MU(4) 

MU(    9)     =  MJ    +    1 

MLM18  )     =  MU(9) 

C0EF(     1)  =    -0. 25D0*ALFA#*2    -    E 

COEF(IO)  =    C0EF(     1) 

COEF(     2)  =    0.5DO*ALFA*(FL    +    0.1D1)    -    Z 

C0EF(12)  =    COEF(    2 ) 

COEF(     3)  =     0.5D0*ALFA*(FN    +    0.1D1)    -    Z 

COEF( 11)  =    CCEF (     3) 

COEFl    4)  =    0.5D0*ALFA*FM 

CCEF{     5)  =    C0EF14) 

C0EFU3)  =  COEFl  4) 

COEF{ 14)  =  CCEF(4) 

COEF (19)  =  -COEF I  4) 

C0EF(2C)  =  -CDEF(  4) 

CGEF(23)  =  -C0EF(4) 

C0EF{24)  =  -COEF (4 ) 

COEF(  C)  =  -FM*(FM  +  FN  +  FL  +  0.1D1) 
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COEF( 15  J  =  COEF(  6) 

COEF(  7)  =  -0.5D0*FL*(FL  +  0.1D1) 

COEF( 17)  =  CGEF{  7) 

COEF{  8  J  =  -0.5DO*FN*(FN  +  0.1D1J 

COEF(lo)  =  CCEF(  8) 

COEF(  9)  =  0.1D1 

COEFU8)  =  COEF(  9) 

COEF(21)  =  FM*FL 

C0EF(26)  =  COEF(21 ) 

COEF(22)  =  FM*FN 

COEF(25J  =  COEF(22) 

IF'(IS.EO.O)  GO  TO  10 

CHANGE  THE  SIGNS  CF  THE  LAST  13  TERMS  FOR  THE  TRIPLET  CASE 
DO  1  1=10,18 

1  CCEF(  IJ  =  -CCEF(  I  ) 
DO  2  1=23,26 

2  CCEF (I )     =    -COEF(  I  ) 
10    LLMAX    =    MAX  I 

NNMAX    =    MA XI    +     1 
HS    =    0.0  DO 
HC    =    O.ODO 
DO    20    K=l,18 

IN    =    NU  { K  J     +    2    +    N N  MAX*  (  L  'i D  A  (  K  )     +    L  L  MAX*  MU  (  K  J     -    1  ) 
HS    =    HS    +    AS(I  N)*COEF  IK] 
HC    =    HC    +    AC( IN)-COEF(K) 
20    CONTINUE 

DC    30    K=19,26 

IN   =    NU(KJ    +   2    +    NNMAX*(LMDA{K)    +    LLMAX*MU(K)    -    1) 

HS    =    HS    +    DS(  I  N)*COEF(  KJ 

HC    =    HC    +    BC( IN)*COEF(K ) 


*i 


RETURN 
END 
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ABSTRACT 


rhc   Harris    expansion    method   is    applied   to   the   elastic    s-wave    scattering    of 
tow    energy    electrons    from    hydrogen    atoms    and   singly    ionized    helium    atoms. 
rhc    trial    wave    functions    are   Hylleraas    functions    of   22,    34,    50   and    70  para- 
meters.      It    is    found    that    reasonably   accurate    values    of   e-H  phase    shifts 
L-an   be   calculated   but    that    e-He+  phase   shifts    are   substantially    less    reliable 
It    is    shown    that    the    Harris    method    gives    an    accurate    depiction    of    the 
location,    but    not    the   width,    of    the    scattering    resonances.      Singlet    and 
triplet    s-wave  phase   shifts    for    e-H   and   e-He+    scattering   are   compared   with 
the   results    of    other    calculations   and  H~   and   He   S   state    energy    levels 
including    the  auto-ionizing    levels    are  presented  and    compared  with    other 
calculations    and   with    experiment.       It    is    tentatively    concluded    that    the 
Harris    method   does    not    work   on    systems    whose    long   range   potential    is    of 
the    Coulomb    type. 
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